(* 
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 = 50G/cm;

RCoil = .846in;
Separation=2.5in;

(* 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 figure 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; (* for copper *)
PackingRatio = 0.9069;
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} (*inches*);

Clear[WireDiameter,Current]

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

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

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

(* Gives good results for Current = 3A, Diameter = 0.0179in or 25AWG *)
Current = 3;
WireDiameter = 0.0254;

{Turns "turns",WireLength, Voltage "V"A/Watts,PowerLoss,CoilWidth"cm"100/m}