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

Update 

* Forgot to include some files in previous commit
parent 300efcdd

matlab/bfield_library_jdl/EFIT_routines/plot_gfile.m

+35 −40
Original line number Diff line number Diff line
@@ -5,17 +5,21 @@ function plot_gfile(g,psi_min,psi_max,npsi,newfig,con_linesty,quiet,plot_psi) 
%

% g can be g struction from readg_g3d, or gfile_name



nRefine = 3;





nRefine = 2;



%% Handle inputs

if ischar(g)

    g = readg_g3d(g);

end

if nargin < 2

    npsi = 50;

end



if nargin < 3 

    psi_min = [];

end

if nargin < 4

    psi_max = [];

end

if nargin < 5

    newfig = 1;

end
@@ -29,26 +33,26 @@ if nargin < 8 
    plot_psi = 'psiN';

end



%% Refine for better contours 

rPlot = g.r;

zPlot = g.z;

for i = 1:nRefine

    [rPlot,zPlot,psiPlot] = refine_psi(rPlot,zPlot,2,g,plot_psi);

end





%% Handle plot_psi switch

switch lower(plot_psi)

    case lower('psiN')

        psiPlot = (g.psirz-g.ssimag)/(g.ssibry-g.ssimag); % Don't need ip_sign here

        sepVal = 1;

        if nargin < 2

        if isempty(psi_min)

            psi_min = 0.6;

        end

        if isempty(psi_max)

            psi_max = 1.1;

        end

        

        cLabel = '\psi_N';

    case lower('psi')

        psiPlot = g.psirz; % Don't need ip_sign here

        sepVal = g.ssibry;

        if nargin < 2

        if isempty(psi_min)

            psi_min = 0.6*(g.ssibry-g.ssimag) + g.ssimag;

        end

        if isempty(psi_max)

            psi_max = 1.1*(g.ssibry-g.ssimag) + g.ssimag;        

        end

        cLabel = '\psi';
@@ -56,29 +60,26 @@ switch lower(plot_psi) 
        error('Bad value for plot_psi')

end







%% Refine for better contours

if nRefine > 1

    rPlot = g.r;

    zPlot = g.z;

    [rPlot,zPlot,psiPlot] = refine_psi(rPlot,zPlot,2^nRefine,g,plot_psi);

end





%%

if newfig == 1

    figure; hold on; box on;

end

% hh=contour(g.r,g.z,psiN_g.',linspace(psi_min,psi_max,npsi),con_linesty,'linewidth',lw);

% % set(hh,'linewidth',2)

% plot(lim_r,lim_z,'k.-')

% plot(bdry_r,bdry_z,'k--')



set(gcf,'color','w');

if isfield(g,'lim')

    plot(g.lim(1,g.lim(1,:)>0),g.lim(2,g.lim(1,:)>0),con_linesty,'linewidth',2)

end

contour(rPlot,zPlot,psiPlot.',linspace(psi_min,psi_max,npsi),'-','linewidth',1);

contour(rPlot,zPlot,psiPlot.',sepVal*[1,1],'k-','linewidth',2);

hc = colorbar;

set(hc,'fontsize',14);

colormap(turbo);

hc.Label.String = cLabel;

hc = colorbar; set(hc,'fontsize',14); colormap(turbo); hc.Label.String = cLabel;

xlabel('R (m)','fontsize',14); ylabel('Z (m)','fontsize',14)

axis equal; axis tight;

ax = axis;

if isfield(g,'gfilename')
@@ -86,7 +87,7 @@ if isfield(g,'gfilename') 
    set(h,'interpreter','none');

end





plot(g.rmaxis,g.zmaxis,'bo')

% Find the xpt(s)

if isfield(g,'bdry')

    xpt_info = find_xpt_jl(g,1,1,1e-8,1);
@@ -94,26 +95,20 @@ if isfield(g,'bdry') 
    xz1 = xpt_info.zx;

    xr2 = xpt_info.rx2;

    xz2 = xpt_info.zx2;



    switch plot_psi

        case lower('psiN')

            contour(rPlot,zPlot,psiPlot.',[1,1]*calc_psiN(g,xr2,xz2),'k-','linewidth',2)



    % contour(g.r,g.z,psiN_g.',[1,1]*2.25,'k-','linewidth',2)



        case lower('psi')

            contour(rPlot,zPlot,psiPlot.',[1,1]*calc_psi(g,xr2,xz2),'k-','linewidth',2)

    end

    plot(xr1,xz1,'bx'); text(xr1+0.01,xz1,'x1','fontsize',8)

    plot(xr2,xz2,'b*'); text(xr2+0.02,xz2,'x2','fontsize',8)

end

% psi_x1 = g.ip_sign*get_psi_bicub(g,xr1,xz1);

% psi_x2 = g.ip_sign*get_psi_bicub(g,xr2,xz2);



% axis

rax = g.rmaxis;

zax = g.zmaxis;

plot(rax,zax,'bo')

xlabel('R [m]','fontsize',14)

ylabel('Z [m]','fontsize',14)

set(gca,'fontsize',14)



end



%%

function [r2,z2,p2] = refine_psi(r,z,fac,g,plot_psi)

r2 = linspace(r(2),r(end-1),length(r)*fac);

z2 = linspace(z(2),z(end-1),length(z)*fac);

matlab/bfield_library_jdl/EFIT_routines/readg_g3d.m

+2 −2
Original line number Diff line number Diff line
@@ -70,8 +70,8 @@ if ~isempty(g.qpsi) % Do not postprocess truncated file 
    g.bicub_coeffs = get_psi_bicub_coeffs(g);

    g.bicub_coeffs_inv = get_psi_bicub_coeffs_inv(g);

    g.fpol_coeffs = polyfit(g.pn,g.fpol,7);

    g.xk = dbsnak(g.mw,g.pn,3);

    g.fpol_coeffs_spline = dbsint(g.mw,g.pn,g.fpol,3,g.xk);

%     g.xk = dbsnak(g.mw,g.pn,3);

%     g.fpol_coeffs_spline = dbsint(g.mw,g.pn,g.fpol,3,g.xk);

    g.filename = filename;

    

    %try to parse filename for shot and time

matlab/bfield_library_jdl/fieldline_following/dl/rk4_core_dl.m

+1 −1
Original line number Diff line number Diff line 
function [yout,ierr] = rk4_core_dl(y,dydx,dx,bfield,nowarn)

if nargin < 6

if nargin < 5

    nowarn = 0;

end

d1 = dx*dydx;

matlab/bjdlTests/calc_B/test_toroidal_flux_g.m

0 → 100644
+60 −0
Original line number Diff line number Diff line 
clearvars;

gfile_name = 'C:\Users\jjl\Dropbox (ORNL)\DIII-D\Haskey\g179477.02280';



g = readg_g3d(gfile_name);





psiThis = calc_psi(g,2.1,0)

% psiThis = -0.1;



% nRefines = -2:7;

% for i = 1:length(nRefines)



nRefine = 4;

% nRefine = nRefines(i)

[rPlot,zPlot,psiPlot] = refine_psi(g.r,g.z,2^nRefine,g);

[r2,z2] = ndgrid(rPlot,zPlot);



in = psiPlot < psiThis;

dA = (rPlot(2)-rPlot(1))*(zPlot(2)-zPlot(1));

b = bfield_geq_bicub(g,r2(in),z2(in));

psi_tor = sum(b.bphi)*dA;





if 0

figure; hold on; box on; set(gcf,'color','w');set(gca,'fontsize',14,'fontweight','bold'); grid on; set(gcf,'color','w');

plot(g.lim(1,g.lim(1,:)>0),g.lim(2,g.lim(1,:)>0),'k-','linewidth',2)

contour(rPlot,zPlot,psiPlot.',psiThis*[1,1],'-','linewidth',1);

contour(rPlot,zPlot,psiPlot.',g.ssibry*[1,1],'k-','linewidth',2);

axis equal; axis tight; ax = axis;

plot(g.rmaxis,g.zmaxis,'bo')

plot(r2,z2,'k.')

plot(r2(in),z2(in),'r.')

end





psiInterp = linspace(g.ssimag,psiThis,100);

psiNInterp = (psiInterp - g.ssimag)/(g.ssibry-g.ssimag);

qInterp = interp1(g.pn,g.qpsi,psiNInterp);

Phi = g.ip_sign*trapz(psiInterp,qInterp)*2*pi





Phi_bry = g.ip_sign*trapz(g.qpsi)*2*pi*(g.ssibry-g.ssimag)/(g.mw-1)



rho = sqrt(Phi/Phi_bry)





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

Reff_bry = sqrt(abs(2*pi*Phi_bry/(pi*g.bcentr)))

rhofromReff = Reff/Reff_bry



% q = dPhi/dPsi



%%

function [r2,z2,p2] = refine_psi(r,z,fac,g)

r2 = linspace(r(2),r(end-1),length(r)*fac);

z2 = linspace(z(2),z(end-1),length(z)*fac);

for i = 1:length(z2)

    p2(:,i) = calc_psi(g,r2,z2(i)*ones(size(r2)));

end

end