Commit 300efcdd authored by Lore, Jeremy's avatar Lore, Jeremy
Browse files

Updates 

*Add simplified mfp routine
*Set NaN for invalid points in Bangle_g
*Add routine for rho from EFIT, rho = sqrt normalized toroidal flux
* Rename old rho routines where rho = R - Rsep
* Cleanup plot_gfile
* Comment out spline routines in gfile
parent 9c1931e1

matlab/bfield_library_jdl/ADAS_routines/calc_mfp.m

0 → 100644
+19 −0
Original line number Diff line number Diff line 
function [mfp,sigma_v] = calc_mfp(te,ne,RateCoeff,indReac,vp)

% Inputs:

%   Te in eV

%   ne in m^-3

%   RateCoeff from read_adas_adf11_file (or equivalent) cm^3/s   (Wcm^3 for plt)

%   indReac: value for last index in RateCoeff.coeff (if scd, ind = charge + 1)

%   vp: test particle velocity (m/s)

% Outputs:

%   mfp in m





ind_oob = find(ne > max(RateCoeff.ne)*1e6);

if ~isempty(ind_oob)

    warning('Some densities are out of bounds in upper direction! Max value used!!')

    ne(ind_oob) = max(RateCoeff.ne)*1e6;

end



sigma_v = 1e-6*10.^(interp2(RateCoeff.te_log10,RateCoeff.ne_log10+6,squeeze(RateCoeff.coeff_log10(:,:,indReac)),log10(te),log10(ne),'linear',NaN));

mfp = vp./(ne.*sigma_v);

matlab/bfield_library_jdl/EFIT_routines/calc_Bangle_g.m

+8 −1
Original line number Diff line number Diff line
@@ -30,6 +30,7 @@ R2 = RZ2(1); 
Z2 = RZ2(2);



DEBUG = 0;

NOWARN = 0;

if DEBUG

    dPlot = 0.01;

    figure; hold on; box on; set(gcf,'color','w');set(gca,'fontsize',14,'fontweight','bold'); grid on;
@@ -64,8 +65,14 @@ for j = 1:length(Reval) 
        t = atan2(vn(2),vn(1));

        plot([Reval(j),Reval(j)+dPlot*cos(t)],[Zeval(j),Zeval(j)+dPlot*sin(t)],'b')        

    end

    b = bfield_geq_bicub(g,Reval(j),Zeval(j),0);

    b = bfield_geq_bicub(g,Reval(j),Zeval(j),NOWARN);

    if isempty(b)

        alpha_deg(j) = NaN;

        beta_deg(j) = NaN;

        continue;

    end

    b = [b.br,b.bz,b.bphi];

       

    b = b./norm(b);

    if DEBUG

        t = atan2(b(2),b(1));

matlab/bfield_library_jdl/EFIT_routines/calc_rho_midplane_from_psiN.m→matlab/bfield_library_jdl/EFIT_routines/calc_R_minus_Rsep_midplane_from_psiN.m

+2 −8
Original line number Diff line number Diff line 
function [rho,R,Rsep] = calc_rho_midplane_from_psiN(g,psiN)



% gfile_name = 'C:\Users\jjl\Dropbox (ORNL)\LDRD 2021\d3d data\g175060.02500';

% 

% g = readg_g3d(gfile_name);



% psiN = [0.99,1.01];

function [RminusRsep,R,Rsep] = calc_R_minus_Rsep_midplane_from_psiN(g,psiN)



nInterp = 100;

rInterp = linspace(g.rmaxis,g.r(end-1),nInterp);
@@ -18,5 +12,5 @@ end 
R = interp1(psiNInterp,rInterp,psiN);

Rsep = interp1(psiNInterp,rInterp,1);



rho = R - Rsep;

RminusRsep = R - Rsep;

matlab/bfield_library_jdl/EFIT_routines/calc_rho_midplane_map.m→matlab/bfield_library_jdl/EFIT_routines/calc_R_minus_Rsep_midplane_map.m

+18 −0
Original line number Diff line number Diff line 
function [rho,R,Rsep] = calc_rho_midplane_map(g,R1,Z1)

% Give R,Z or just R = psiN

% rho = R - Rsep



if nargin == 2

    psiN1 = R1;

end

function [RminusRsep,R,Rsep] = calc_R_minus_Rsep_midplane_map(g,R1,Z1)

% Give R,Z 

% RminusRsep = R - Rsep



Ninterp = 100;

Rinterp = linspace(g.rmaxis,g.r(end-1),Ninterp);

Zinterp = g.zmaxis*ones(1,Ninterp);

[psiN_mid,psi_mid] = calc_psiN(g,Rinterp,Zinterp);



% plot_gfile(g)

% plot(Rinterp,Zinterp,'cx')

% plot(R1,Z1,'cx','markersize',12,'linewidth',3)





if nargin == 3

[psiN_mid,psi_mid] = calc_psiN(g,Rinterp,Zinterp);

[psiN1,psi1] = calc_psiN(g,R1,Z1);

end



R_midplane_map = interp1(psiN_mid,Rinterp,psiN1);   



Rsep = interp1(psiN_mid,Rinterp,1);   

rho_midplane_map = R_midplane_map - Rsep;



rho = rho_midplane_map;

RminusRsep = rho_midplane_map;

R = R_midplane_map;
 


% figure; hold on;

% plot(Rinterp,psiN_mid)

% plot(R_midplane_map,psiN1,'ro')

matlab/bfield_library_jdl/EFIT_routines/calc_rho.m

0 → 100644
+23 −0
Original line number Diff line number Diff line 
function [rho,Phi,Reff,Reff_bry] = calc_rho(g,R1,Z1)

% [rho] = calc_rho(g,R1,Z1)

% Rho is sqrt toroidal flux

% Phi is toroidal flux

% Reff is R from Bo*pi*R^2 = 2*pi*Phi



% See test_toroidal_flux 



%% Map to psiN first

psiN = calc_psiN(g,R1,Z1);

if psiN >= 1

    error('rho undefined outside of separatrix. PsiN here = %f\n',psiN)

end



% Integrate q = Phi'/Psi' and interpolate at psiN

PhiEval = g.ip_sign*cumtrapz(g.qpsi)*2*pi*(g.ssibry-g.ssimag)/(g.mw-1);

Phi = interp1(g.pn,PhiEval,psiN);

rho = sqrt(Phi/PhiEval(end));

Reff = sqrt(abs(2*pi*Phi/(pi*g.bcentr)));

Reff_bry = sqrt(abs(2*pi*PhiEval(end)/(pi*g.bcentr)));