Commit ad0373fb authored by Lore, Jeremy's avatar Lore, Jeremy
Browse files

Major update to gfile interpolation routines 

*Behavior may change!
*Results now uniform in region of applicability, nans returned when in boundary cells or outside
*add new calc_bicub_interpolation_inds
*Remove old bicub routines, replace with inverse bicub
*Calls to calc_psi, calc_psiN, bfield_geq_bicub replaced with updated routines. These will give the same result, but points off grid are now nan instead of routines returning empty.
*get_psi_bicub removed, use calc_psi instead!
*bugfix in xpoint evaluation (does not affect output, just allows to run for different resolutions)
*cleanup plot_gfile, some old routines
parent 91c8d656

matlab/bfield_library_jdl/EFIT_routines/calc_bicub_interpolation_inds.m

0 → 100644
+31 −0
Original line number Diff line number Diff line 
function [ir,iz,index,ierr] = calc_bicub_interpolation_inds(g,R,Z,nowarn)

% Find the position in the psirz grid and return indices for interpolation



if iscolumn(R)

    R = R';

    Z = Z';

end



if any(isnan(R)) || any(isnan(Z))

    error('R or Z is nan')

end



% Add a small offset to avoid error with first grid point

ir = floor( (R - g.r(1) + 1e-10)/g.dR ) + 1;  

iz = floor( (Z - g.z(1) + 1e-10)/g.dZ ) + 1;



% check for points off grid, no derivatives on boundary cells

ierr = false(size(R));

ierr(ir < 2 | ir > g.mw - 1) = true;

ierr(iz < 2 | iz > g.mh - 1) = true;



ir(ierr) = 1;

iz(ierr) = 1;



index = ir + g.mw*(iz-1);



if ~nowarn

    for i = find(ierr)

        warning('Point off grid: [R,Z] = [%.2f,%.2f], [Rmin,Zmin] = [%.2f,%.2f], [Rmax,Zmax] = [%.2f,%.2f]\n',R(ierr(i)),Z(ierr(i)),g.r(2),g.r(g.mw-1),g.z(2),g.z(g.mh-1))

    end

end
 
 No newline at end of file

matlab/bfield_library_jdl/EFIT_routines/calc_psi.m

+72 −7
Original line number Diff line number Diff line 
function [psi] = calc_psi(g,R1,Z1,quiet)

%[psi] = calc_psi(g,R1,Z1,quiet)

% R1, Z1 can be vectors

% quiet optional  

%--- Just applies g.ip_sign to get_psi_bicub!

function [psi,ierr,dpsidr,dpsidz] = calc_psi(g,R,Z,nowarn)



if nargin < 4

    quiet = 0;

    nowarn = 0;

end



npts = length(R);

R=reshape(R,1,npts);

Z=reshape(Z,1,npts);



[ir,iz,index,ierr] = calc_bicub_interpolation_inds(g,R,Z,nowarn);



dr = (R - g.r(ir))/g.dR;

dz = (Z - g.z(iz))/g.dZ;



c = g.bicub_coeffs_inv;



c00 = c.c00(index);

c10 = c.c10(index);

c20 = c.c20(index);

c30 = c.c30(index);



c01 = c.c01(index);

c11 = c.c11(index);

c21 = c.c21(index);

c31 = c.c31(index);



c02 = c.c02(index);

c12 = c.c12(index);

c22 = c.c22(index);

c32 = c.c32(index);



c03 = c.c03(index);

c13 = c.c13(index);

c23 = c.c23(index);

c33 = c.c33(index);



drr = dr.*dr;

drrr = drr.*dr;

dzz = dz.*dz;

dzzz = dzz.*dz;



drz   = dr.*dz;

drzz  = dr.*dzz;

drzzz = dr.*dzzz;

drrz  = drr.*dz;

drrzz = drr.*dzz;

drrzzz = drr.*dzzz;

drrrz  = drrr.*dz;

drrrzz = drrr.*dzz;



psi = g.ip_sign*(...

    c00       + c10.*dr    + c20.*drr    + c30.*drrr   + ...

    c01.*dz   + c11.*drz   + c21.*drrz   + c31.*drrrz  + ...

    c02.*dzz  + c12.*drzz  + c22.*drrzz  + c32.*drrrzz + ...

    c03.*dzzz + c13.*drzzz + c23.*drrzzz + c33.*drrr.*dzzz);



psi(ierr) = nan;



if nargin > 2

    dpsidr = g.ip_sign*(...

        c10       + 2*c20.*dr    + 3*c30.*drr   + ...

        c11.*dz   + 2*c21.*drz   + 3*c31.*drrz  + ...

        c12.*dzz  + 2*c22.*drzz  + 3*c32.*drrzz + ...

        c13.*dzzz + 2*c23.*drzzz + 3*c33.*drrzzz)./g.dR;



    dpsidz = g.ip_sign*(...

          c01      +   c11.*dr   +   c21.*drr   +   c31.*drrr   + ...

        2*c02.*dz  + 2*c12.*drz  + 2*c22.*drrz  + 2*c32.*drrrz  + ...

        3*c03.*dzz + 3*c13.*drzz + 3*c23.*drrzz + 3*c33.*drrrzz)./g.dZ;



    dpsidr(ierr) = nan;

    dpsidz(ierr) = nan;

end
 
 No newline at end of file
 
psi = g.ip_sign*get_psi_bicub(g,R1,Z1,quiet);
 
 No newline at end of file

matlab/bfield_library_jdl/EFIT_routines/calc_psiN.m

+6 −5
Original line number Diff line number Diff line 
function [psiN,psi] = calc_psiN(g,R1,Z1,quiet)

%[psiN,psi] = calc_psiN(g,R1,Z1,quiet)

% R1, Z1 can be vectors

function [psiN,psi,ierr] = calc_psiN(g,R,Z,quiet)

%[psiN,psi] = calc_psiN(g,R,Z,quiet)

% R, Z can be vectors

% quiet optional



if nargin < 4

    quiet = 0;

end

psi = g.ip_sign*get_psi_bicub(g,R1,Z1,quiet);

psiN = (psi-g.ssimag)/(g.ssibry-g.ssimag);
 
 No newline at end of file
 
[psi,ierr]= calc_psi(g,R,Z,quiet);         % this applies g.ip_sign

psiN = (g.ip_sign*psi-g.ssimag)/(g.ssibry-g.ssimag); % so have to apply it again because ssi quantities will be flipped (if ip_sign = -1)
 
 No newline at end of file

matlab/bfield_library_jdl/EFIT_routines/find_xpt_jl.m

+4 −3
Original line number Diff line number Diff line
@@ -93,8 +93,8 @@ if refine == 1 
  zg = zeros(n1,n1);

  niter = 1;

  while err > tol && niter < niter_max

    rt = linspace(max([Rmin_eval,rx-der]),min([Rmax_eval,rx+der,n1]));

    zt = linspace(max([Zmin_eval,zx-dez]),min([Zmax_eval,zx+dez,n1]));

    rt = linspace(max([Rmin_eval,rx-der]),min([Rmax_eval,rx+der]),n1);

    zt = linspace(max([Zmin_eval,zx-dez]),min([Zmax_eval,zx+dez]),n1);

    for i = 1:n1 

      ztmp = repmat(zt(i),1,n1);

      b = bfield_geq_bicub(g,rt,ztmp);
@@ -115,6 +115,7 @@ if refine == 1 
    zx = zg(ix,jx);



    niter = niter + 1;    

    

    if niter >= niter_max

        warning('niter_max exceeded for 1st x-point')

    end

matlab/bfield_library_jdl/EFIT_routines/get_bicub_mat.m

deleted100644 → 0
+0 −36
Original line number Diff line number Diff line 
function bicub_mat = get_bicub_mat()



bicub_mat = zeros(16,16);

% ;;;;;; function values at corners

bicub_mat(1,1) = 1.0;

bicub_mat(2,[1,2,3,4]) = [1.,1.,1.,1.];

bicub_mat(3,[1,5,9,13]) = [1.,1.,1.,1.];

bicub_mat(4,:) = 1.;



% ;;;;;;; 1st derivatives at corners: x direction

bicub_mat(5,2)=1.0;

bicub_mat(6,[2,3,4]) = [1.0,2.0,3.0];

bicub_mat(7,[2,6,10,14]) = [1.0,1.0,1.0,1.0];

bicub_mat(8,[2,3,4]) = [1.,2.,3.];

bicub_mat(8,[6,7,8]) = [1.,2.,3.];

bicub_mat(8,[10,11,12]) = [1.,2.,3.];

bicub_mat(8,[14,15,16]) = [1.,2.,3.];



% ;;;;;;; 1st derivatives at corners: y direction

bicub_mat(9,5) = 1.0;

bicub_mat(10,5:8) = [1.,1.,1.,1.];

bicub_mat(11,[5,9,13]) = [1.,2.,3.];

bicub_mat(12,[5,9,13]) = [1.,2.,3.];

bicub_mat(12,[6,10,14]) = [1.,2.,3.];

bicub_mat(12,[7,11,15]) = [1.,2.,3.];

bicub_mat(12,[8,12,16]) = [1.,2.,3.];



% ;;;;;;; cross derivatives at corners

bicub_mat(13,6) = 1.;

bicub_mat(14,[6,7,8]) = [1.,2.,3.];

bicub_mat(15,[6,10,14]) = [1.,2.,3.];

bicub_mat(16,[6,10,14]) = [1.,2.,3.];

bicub_mat(16,[7,11,15]) = [2.,4.,6.];

bicub_mat(16,[8,12,16]) =[3.,6.,9.];



% bicub_mat = bicub_mat.';
 
 No newline at end of file