clc
clear all
%Setup up: This program calculates the inductance, resistance, and mutual
%inductance between two multi-layered coils. The coils are composed of
%wire with circular cross-sections (as opposed to wire of rectangular cross
%section that you would find on a PCB). 
%%%% Initialize, clear out workspace variables.
%%%%%%%%%%%%%%%USER INPUT VALUES
f=7E5; %frequency of operationc [Hz]
cond=5.8E7; %conductivity of wires
rw = 2.5527e-04; %wire radius [meters]
             %AWG 20 -> rw = 4.0640e-04
             %AWG 21 -> rw = 3.6195e-04
             %AWG 22 -> rw = 3.2258e-04
             %AWG 24 -> rw = 2.5527e-04
Np=10;%number of turns (primary coil on generator side)
Ns=10; %number of turns (secondary coil on load side)
gzp=3*rw; %Winding pitch of primary coil (i.e. center to center spacing between windings of primary coil... should be larger than 2*rw... 3*rw is a good estimate for magnet wire) [meters] 
gzs=3*rw;%Winding pitch of secondary coil (i.e. center to center spacing between windings of secondary coil.. should be larger than 2*rw... 3*rw is a good estimate for magnet wire) [meters] 
Rp=.161925; %Radius of the coil (primary) [meters]
Rs=.161925; %Radius of the coil (secondary) [meters]
% Radius to external surface of 12 inch Sch.40 PVC pipe
dist=0.1; %distance seperating primary and secodary spirals, i.e.
%PRIMARY SPRIAL  <----dist-----> SECONDARY SPIRAL [meters]
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%CODE BEGIN!
zp=[-(Np-1)*gzp:gzp:0]; %z-positions of loops in primary coil
zs=[dist:gzp:dist+(Ns-1)*gzp]; %z-positions of loops in primary coil
%FIND SELF INDUCTANCE OF PRIMARY LOOP:
%calculate mutual inductance between every loop in the primary spriral 
for n=1:length(zp)
    for m=1:length(zp)
        z_distance_between=abs(zp(n)-zp(m)); %calculate distance between loops
        Lp_circ(n,m)=coaxial_circle(Rp,Rp,rw,z_distance_between); %calculate inductance between loops
    end
end
Primary_inductance_circular_loops = sum(sum(Lp_circ)); %Sum the contributions
                   %from all loops to find total self-inductance of primary
%FIND SELF INDUCTANCE OF SECONDARY LOOP:
%calculate mutual inductance between every loop in the primary spriral 
for n=1:length(zs)
    for m=1:length(zs)
        z_distance_between=abs(zs(n)-zs(m)); %calculate distance between loops
        Ls_circ(n,m)=coaxial_circle(Rs,Rs,rw,z_distance_between); %calculate inductance between loops
    end
end
Secondary_inductance_circular_loops = sum(sum(Ls_circ)); %Sum the contributions
                   %from all loops to find total self-inductance of primary
%calculate mutual inductance
for n=1:length(zp)
    for m=1:length(zs)
        z_distance_between=abs(zp(n))+abs(zs(m)); %distance between each loop
        M_circ(n,m)=coaxial_circle(Rp,Rs,rw,z_distance_between); %calculate inductance between loops
    end
end
Mutual_inductance_circular_loops = sum(sum(M_circ)); %Sum the contributions 
                 %from all loops to find total mutual inductance
%%%%%%%Find resistances of each loop (circular loops)
skin_depth=sqrt(1/cond/pi/4/pi/10^-7/f);%skin depth
%(fron "A simple formula for calculating the
%frequency-dependent resitance of a round wire"... in this reference equation (15) is correct... equation (16) is off by factor of 2!)
R_skin=1/cond/skin_depth;
resistance_primary=1/cond/pi/skin_depth/(1-exp(-rw/skin_depth))/(2*rw-skin_depth*(1-exp(-rw/skin_depth)))*2*pi*Rp*Np;
resistance_secondary=1/cond/pi/skin_depth/(1-exp(-rw/skin_depth))/(2*rw-skin_depth*(1-exp(-rw/skin_depth)))*2*pi*Rs*Ns;            
%calculate C
Cp=1/Primary_inductance_circular_loops/(2*pi*f)^2;
Cs=1/Secondary_inductance_circular_loops/(2*pi*f)^2;
%calculate U (from Witricity paper)
U=2*pi*f*Mutual_inductance_circular_loops/sqrt(resistance_secondary*resistance_primary);
%calculate ideal generator impedance
Rg=resistance_primary*sqrt(1+U^2);
%calculate ideal load impedance
RL=resistance_secondary*sqrt(1+U^2);
%calculate maximum efficiency
eta_max=U^2/(1+sqrt(1+U^2))^2;
%calculate critical coefficient of coupling (for lossless case with same
%generator and load resistance values Rg and RL). 
Q1=(f*2*pi*Primary_inductance_circular_loops)/Rg;
Q2=(2*pi*f*Secondary_inductance_circular_loops)/RL;
Kcrit=1/(sqrt(Q1*Q2));
disp('Results:')
disp(' ')
disp(['The primary inductance is ', num2str(Primary_inductance_circular_loops), ' [H].'])
disp(['The secondary inductance is ', num2str(Secondary_inductance_circular_loops), ' [H].'])
disp(['The mutual inductance between secondary and primary loops is ', num2str(Mutual_inductance_circular_loops), ' [H].'])
disp(['The coupling coefficient k is ', num2str(Mutual_inductance_circular_loops/sqrt(Primary_inductance_circular_loops*Secondary_inductance_circular_loops))])
disp(['The resistance of the primary loop is ',num2str(resistance_primary), ' [Ohms].'])
disp(['The resistance of the secondary loop is ',num2str(resistance_secondary), ' [Ohms].'])
disp(['The required capactiance in the primary loop for resonance is ',num2str(Cp), ' [F].'])
disp(['The required capactiance in the secondary loop for resonance is ',num2str(Cs), ' [F].'])
disp(['The ideal generator resistance for maximum attainable power transmission* is ',num2str(Rg), ' [Ohms].'])
disp(['The ideal load resistance for maximum attainable power transmission* is ',num2str(RL), ' [Ohms].'])
disp(['The maximum attainable power transmission* is ',num2str(eta_max*100), ' [%].'])
disp(['*Here, power transmission is defined as the ratio: (Power delivered to load)/(Maximum power available from generator)'])
%disp(['The critical coupling coefficient is ',num2str(Kcrit)])
%