## Copyright (C) 2009 Indrek Mandre <indrek(at)mare.ee>
## For more information please see http://www.mare.ee/indrek/octave/

% Parameters available:

global Rp = 0.15;                         % physical coil radius
global Sp = 0.08;                         % physical coil spacing
global Ip = 100000;                       % physical coil current

global A = 0.15;                          % desired wiffleball radius

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
inv_f = (A^2)/(Rp^2 + (Rp + Sp/sqrt(2))^2);  % spherical inversion factor
global Ri = Rp * inv_f ;                  % image coil radius
global Si = Sp * inv_f ;                  % image coil spacing
global Ii = -Ip / sqrt(inv_f);            % image coil current

global rr = Rp + Sp / sqrt(2);
global clen = 2 * pi * Rp;                % coil length

global dowb = false;                      % add WB or not into force_z()

% calculate field from a polywell, with the top coil missing
function [B] = polywell_5coils(radius, spacing, current, pos)
  rr = radius + spacing / sqrt(2);
  B = loop_bfield3d([-rr,   0,    0], [ 1,   0,    0], radius, current, pos) + \
      loop_bfield3d([ rr,   0,    0], [-1,   0,    0], radius, current, pos) + \
      loop_bfield3d([  0, -rr,    0], [ 0,   1,    0], radius, current, pos) + \
      loop_bfield3d([  0,  rr,    0], [ 0,  -1,    0], radius, current, pos) + \
      loop_bfield3d([  0,   0,  -rr], [ 0,   0,    1], radius, current, pos);
endfunction

function y = force_z(x)
  global rr clen Rp Sp Ip Ri Si Ii dowb
  angle = (x / clen) * 2 * pi;
  pos = [ Rp * cos(angle), Rp * sin(angle), rr];
  L = -[ -sin(angle), cos(angle), 0 ];
  B = polywell_5coils(Rp, Sp, Ip, pos);
  % and add the image coils
  if dowb
    B = B .+ polywell_bfield (Ri, Si, Ii, pos);
  endif
  F = Ip .* cross(L,B);
  y = F(3); % take the z component
endfunction

% Calculate the force on a single coil
function f = calc_force(radius, spacing, current, wbradius = 0)
  global rr clen Rp Sp Ip Ri Si Ii dowb
  Rp = radius;
  Sp = spacing;
  Ip = current;

  rr = Rp + Sp / sqrt(2);
  clen = 2 * pi * Rp;

  dowb = (wbradius != 0);
  if dowb
    A = wbradius;
    inv_f = (A^2)/(Rp^2 + (Rp + Sp/sqrt(2))^2);
    Ri = Rp * inv_f;
    Si = Sp * inv_f;
    Ii = -Ip / sqrt(inv_f);
  endif

  f = quad (@force_z, 0, clen);
endfunction

% Calculating the WB pressure for R=0.15, S=0.08 with various currents
% and wiffleball radiuses.
data = zeros(size(0.1:0.0125:0.2,2),size(1:5,2) + 1);
f1amp = calc_force(0.15, 0.08, 1, 0);

for r=1:16
  wbr = 0.05 + (r - 1) * (0.2 - 0.05) / 15;
  data(r,1) = wbr;
  wbf1amp = calc_force (0.15, 0.08, 1, wbr);
  for i=1:5
    current = i * 100000;
    wbf = wbf1amp * current ** 2;
    f = f1amp * current ** 2;
    df = 6 * (wbf - f);
    sa = 4 * pi * wbr ** 2;
    p = df / sa;
    data(r,i+1) = p;
    printf ("I=%f r=%.4f P=%.2f Pa\n", current, wbr, p);
  endfor
endfor

data
save ("-ascii", "polywell_force_wbp.dat", "data");

% Calculating the forces on coils dependance on current
% for various R/S ratios.
force = zeros(5,1+6);

for r=1:6
  ratio = 1.375 + (r-1) * 0.250;
  ratio
  f1amp = calc_force (0.15, 0.15 / ratio, 1);
  for i=1:5
    Ip = i * i * 50000;
    force(i,1) = Ip / 1000;
    force(i,1 + r) = f1amp*Ip**2 / 1000;
  endfor
endfor

force
save ("-ascii", "polywell_force_force.dat", "force");

RADIUS = 0.15;
SPACING = 0.08;
wbf = zeros (14,6);

for r = 1:14
  ratio = (r-1) * (1 / 10.0);
  f1amp = calc_force (RADIUS, SPACING, 1);
  wbf1amp = calc_force (RADIUS, SPACING, 1, ratio * RADIUS);
  wbf(r,1) = ratio;
  for i=1:5
    Ip = i * 200000;
    wbf(r,i + 1) = (wbf1amp * Ip ** 2 - f1amp * Ip ** 2) / 1000;
    printf ("RADIUS=%.2f, SPACING=%.2f, A/R=%.3f, CURRENT=%.0f, WBF=%.2f\n", RADIUS, SPACING, ratio, Ip, wbf(r,i + 1));
  endfor
endfor

wbf
save ("-ascii", "polywell_force_wbf.dat", "wbf");