(* 
This notebook describes the magnetic field of an Anti-helmholtz coil for the Strontium 6in Spherical Octagon along the axis of the coil. The formula was derived using the Biot-Savart Law. 

Sorrentino gives the magnetic field gradient near the center of their MOT as 50 Gauss/cm.

However, Dan says that we'll only need ~15 Gauss/cm.
*)

(* Definitions *)
in = 2.54cm;
cm = 0.01m;
T = 10000G;

(* Constants *)
μ0 = 4Pi 10^-7 T m/A;

(* Variables *)
RequiredGradient = 15G/cm;

RCoil = 7.2/2in;
Separation=4.430in;

(* Calculate the B field along the axis of symmetry. This field is in the z direction. The quantity "CurrentTurns" is what this calculation aims to determine, which depends only on geometery and the RequiredGradient, which we will then use to find wire diameter, turns, etc *)

Clear[z,CurrentTurns]

B =μ0 CurrentTurns (RCoil/m)^2 /(2 m)(((RCoil/m)^2+(z-Separation/(2 m))^2)^(-3/2)-((RCoil/m)^2+(z+Separation/(2 m))^2)^(-3/2));

BGradient = D[B,z]/m;

z=0;
CurrentTurns = (CurrentTurns/.Solve[A BGradient==  RequiredGradient ,CurrentTurns][[1]])A;
Clear[z];

{CurrentTurns"CurrentTurns"/A,

Plot[B/G,{z,-2RCoil/m,2RCoil/m},AxesLabel->{"x (m)","B (G)"},ImageSize->Small],

Plot[BGradient m/G/100,{z,-2RCoil/m,2RCoil/m},AxesLabel->{"x (m)","dB/dz (G/cm)"},ImageSize->Small]}

(* Now that we have a value for CurrentTurns, we may find wire parameters that will produce it with reasonable Ohmic losses and size *)

(* Definitions and Constants *)
Ω = Watts/A^2 ;
(* ρ = 1.68 10^-8 Ω m;*) (* Theoretical resistivity of copper *)
ρ = 1.8558 10^-8 Ω m ;        (* Advertized by MWS *)
PackingRatio = 0.9069;   (* Area of hcp circles over total area *)
InsulationThickness = (.0335-.0320)in;
AWG = {.2893(*1*),.2576,.2294,.2043,.1819,.162,.1443,.1285,.1144,.1019 (*10*),.0907,.0808,.072,.0641,.0571,.0508,.0453,.0403,.0359,.032(*20*),.0285,.0254,.0226,.0201,.0179(*25*),.0159,.0142,.0126,.0113,.01 (*30*)} (*inches*);

Clear[WireDiameter,Current]

(* Calculate the ohmic losses *)
Turns = CurrentTurns/(Current A);
WireLength = 2 RCoil Pi Turns;
WireArea = Pi (WireDiameter in)^2/4;
InsulatedArea = Pi(WireDiameter in+InsulationThickness)^2/4;
Resistance = ρ WireLength/WireArea;
Voltage = Current A Resistance; (* Voltage is independent of Current here! *)
PowerLoss = Current A Voltage;

(* And the dimensions of the coil, assuming square cross section with hexagonal closest packing *)
CoilArea = Turns InsulatedArea/PackingRatio;
CoilWidth =Simplify[CoilArea^(1/2),m>0];

Animate[
Evaluate[
{{Current "A",WireDiameter "in"},{Turns "turns",WireLength, Voltage "V"A/Watts,PowerLoss,CoilArea"cm^2"10000/m^2,CoilWidth"cm"100/m}}
],
{Current,1,10},{WireDiameter,AWG},AnimationRunning->False]

(* Coil area must be less than 2.322 cm^2. Current must be less than 10A. Power must be less than 100Watts. It would be nice to minimize number of turns *)

(* Gives good results for Current = 3A, Diameter = 0.0285in or 21AWG *)

Current = 4; 
WireDiameter = 0.032;


{Turns "turns",WireLength, Voltage "V"A/Watts,PowerLoss,CoilArea"cm^2"10000/m^2,CoilArea/(.5*2.54^2)"in height"10000/m^2}