% conpute effective lengths for hitsi
% first get boundary data
% bdata(:,1) = r [m]
% bdata(:,2) = z [m]
% bdata(:,3) = Bp [A/m]
% bdata(:,4) = l (arclength) [m]
% bdata(:,5) = current function [A]
% I_tor = 2 * bdata(np,5)
clear
bb=dlmread('bb.txt');
r0=max(bb(:,1))/2;
bt=atan2(bb(:,2),r0-bb(:,1));
[a,ia,ja]=unique(bt);
bs=bb(ia,1:3);
top=bs(:,2)>=0;
np=sum(top)
bdata=zeros(np,5);
bdata(:,1:3)=bs(top,:);
for i=2:np
    dl=sqrt(sum((bdata(i,1:2)-bdata(i-1,1:2)).^2));
    ba=(bdata(i,3)+bdata(i-1,3))/2;
    bdata(i,4)=bdata(i-1,4)+dl;
    bdata(i,5)=bdata(i-1,5)+dl*ba;
end
% normalize to unit current & MKSA units
mu0=4.e-7*pi;
tcur=2*bdata(np,5);
bdata(:,3)=bdata(:,3)*mu0/tcur;
bdata(:,5)=bdata(:,5)/tcur;
plot(bdata(:,4),bdata(:,3))

% hitsi flux loops
fl.r=[.08274;...
      .11866;...
      .15458;...
      .19050;...
      .22643;...
      .41927;...
      .44235;...
      .46542;...
      .48850;...
      .51157];
fl.n=['FS05 & FS06';...
      'FS04 & FS07';...
      'FS03 & FS08';...
      'FS02 & FS09';...
      'FS01 & FS10';...
      'FL01 & FL10';...
      'FL02 & FL09';...
      'FL03 & FL08';...
      'FL04 & FL07';...
      'FL05 & FL06'];
nlp=length(fl.r)
for i=1:nlp
    fl.l(i,1)=fzero(@(l)(interp1(bdata(:,4),...
    bdata(:,1),l)-fl.r(i,1)),[0,bdata(np,4)]);
end
fl.z=interp1(bdata(:,4),bdata(:,2),fl.l);
fl.b=interp1(bdata(:,4),bdata(:,3),fl.l);
brk=zeros(nlp+1,1);
brk(nlp+1,1)=bdata(np,4);
for i=2:nlp
    brk(i,1)=(fl.l(i,1)+fl.l(i-1,1))/2;
end
cur=interp1(bdata(:,4),bdata(:,5),brk);
fl.c=cur(2:nlp+1,1)-cur(1:nlp,1);
fl.d=mu0*fl.c./fl.b;
for i=1:nlp
    disp(sprintf([fl.n(i,:),'  %f  %f  %e'],...
        [fl.l(i,1),fl.d(i,1),fl.b(i,1)]))
end

% hitsi surface probes (minus gap)
sp.r=[.10060;...
      .13652;...
      .17254;...
      .20837;...
      .43081;...
      .45388;...
      .47696;...
      .50003];
%      .54928];
sp.n=['S04P & S05P';...
      'S03P & S06P';...
      'S02P & S07P';...
      'S01P & S08P';...
      'L01P & L10P';...
      'L02P & L09P';...
      'L03P & L08P';...
      'L04P & L07P'];
%      'L05P & L06P'];
npr=length(sp.r)
for i=1:npr
    sp.l(i,1)=fzero(@(l)(interp1(bdata(:,4),...
    bdata(:,1),l)-sp.r(i,1)),[0,bdata(np,4)]);
end
sp.z=interp1(bdata(:,4),bdata(:,2),sp.l);
sp.b=interp1(bdata(:,4),bdata(:,3),sp.l);
brk=zeros(npr+1,1);
brk(npr+1,1)=bdata(np,4);
for i=2:npr
    brk(i,1)=(sp.l(i,1)+sp.l(i-1,1))/2;
end
cur=interp1(bdata(:,4),bdata(:,5),brk);
sp.c=cur(2:npr+1,1)-cur(1:npr,1);
sp.d=mu0*sp.c./sp.b;
for i=1:npr
    disp(sprintf([sp.n(i,:),'  %f  %f  %e'],...
        [sp.l(i,1),sp.d(i,1),sp.b(i,1)]))
end