equilibrium_jd

PURPOSE ^

SYNOPSIS ^

function [varargout] = equilibrium_jd(equil,radialDKE,display,method,mfactor,axs,font)

DESCRIPTION ^

   This function interpolates the equilibrium data on the radial grid.

   It first performs a very precise FFT-based interpolation in theta in
   order to obtain precise future interpolation and integration
   calculations in theta.

   It also calculates several quantities linked to the equilibrium, such
   as the trapped/passing boundary in momentum space, the plasma minor
   radius on LFS midplane, and the effective charge. 

   NOTE: The equilibrium data in EQUIL_tokname must follow these
   conventions:

       - R Z BR BZ BPHI must be given on a (psi,theta) grid [npsi,ntheta].
       - ppsi(1) = 0 is the magnetic axis.
       - ppsi(npsi) is the poloidal flux on the last closed flux-surface.
           It can be either positive or negative depending on the sign of the
           current with respect to the toroidal angle phi.
       - theta(1) = 0 and theta(ntheta) = 2pi is required.
       - The psi can be non-uniform. The theta grid MUST be uniform
       - other equilibrium parameters (temperatures, densities) are
           assumed to be uniform within a flux-surface, because of the fast parallel
           motion. They are therefore given on the psi grid.

   INPUTS:

       - tokamak: tokamak scenario selected
       - xpsin: positions for FP calculations on the normalized psi grid
           (NAN: ppsi grid)
       - xpsin_S_dke: radial flux grid used for the calculation of FS volume and poloidal surface 
       - display_mode: (0-1) text only (2) graphics
       - method: interpolation method (linear,pchip,spline,...) 
           (default: linear)
       - mfactor: factor for super theta grid (default: 50)
       - axs: list of axes to update
       - font : structure containing font parameters for graphs

   OUTPUTS:

       - theta: poloidal angle grid [1,ntheta]
       - xpsin: normalized poloidal flux grid [1,nr]
       - psia_apRp: Edge normalized cylindrical poloidal flux (Wb) [1,1]
       - Xx: major radius [nr,ntheta]
       - Xy: elevation [nr,ntheta]
       - XBx: B-field component in R direction [nr,ntheta]
       - XBy: B-field component in Z direction [nr,ntheta]
       - XBPHI: B-field component in PHI direction [nr,ntheta]
       - xB0: minimum magnetic field on a flux-surface [1,nr]
       - xR0: corresponding major radius [1,nr]
       - xBT0: corresponding toroidal field [1,nr]
       - xBp0: corresponding poloidal field [1,nr]
       - xrho: normalized radius on the horizontal LFS midplane [1,nr]
       - xTe: electron temperature [1,nr]
       - xne: electron density [1,nr]
       - xzTi: ion temperatures [ni,nr]
       - xzni: ion densities [ni,nr]
       - zZi: ion charges [1,ni]
       - zmi: ion masses [1,ni]
       - xZeff: effective charge [1,nr]
       - Rp: major radius on axis [1,1]
       - Zp: elevation on axis [1,1]
       - Bax: magnetic field on axis [1,1]
       - ap: minor radius on LFS midplane [1,1]
       - Ip: Toroidal plasma current [1,1]
       - Ncoils: Number of toroidal magnetic field coils [1,1]
       - BeLi: Beta + li/2 [1,1]
       - xmhu0T2: pitch-angle (squared) at the trapped/passing boundary [1,nr]
       - xdV_dke: incremental FS volume at dke positions
       - xdA_dke: incremental FS poloidal surface at dke positions
       - HXR: Fast electron bremsstrahlung

   By J. Decker RLE-MIT <jodecker@alum.mit.edu> and Y. Peysson DRFC-CEA <yves.peysson@cea.fr>

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function [varargout] = equilibrium_jd(equil,radialDKE,display,method,mfactor,axs,font)
0002 %
0003 %   This function interpolates the equilibrium data on the radial grid.
0004 %
0005 %   It first performs a very precise FFT-based interpolation in theta in
0006 %   order to obtain precise future interpolation and integration
0007 %   calculations in theta.
0008 %
0009 %   It also calculates several quantities linked to the equilibrium, such
0010 %   as the trapped/passing boundary in momentum space, the plasma minor
0011 %   radius on LFS midplane, and the effective charge.
0012 %
0013 %   NOTE: The equilibrium data in EQUIL_tokname must follow these
0014 %   conventions:
0015 %
0016 %       - R Z BR BZ BPHI must be given on a (psi,theta) grid [npsi,ntheta].
0017 %       - ppsi(1) = 0 is the magnetic axis.
0018 %       - ppsi(npsi) is the poloidal flux on the last closed flux-surface.
0019 %           It can be either positive or negative depending on the sign of the
0020 %           current with respect to the toroidal angle phi.
0021 %       - theta(1) = 0 and theta(ntheta) = 2pi is required.
0022 %       - The psi can be non-uniform. The theta grid MUST be uniform
0023 %       - other equilibrium parameters (temperatures, densities) are
0024 %           assumed to be uniform within a flux-surface, because of the fast parallel
0025 %           motion. They are therefore given on the psi grid.
0026 %
0027 %   INPUTS:
0028 %
0029 %       - tokamak: tokamak scenario selected
0030 %       - xpsin: positions for FP calculations on the normalized psi grid
0031 %           (NAN: ppsi grid)
0032 %       - xpsin_S_dke: radial flux grid used for the calculation of FS volume and poloidal surface
0033 %       - display_mode: (0-1) text only (2) graphics
0034 %       - method: interpolation method (linear,pchip,spline,...)
0035 %           (default: linear)
0036 %       - mfactor: factor for super theta grid (default: 50)
0037 %       - axs: list of axes to update
0038 %       - font : structure containing font parameters for graphs
0039 %
0040 %   OUTPUTS:
0041 %
0042 %       - theta: poloidal angle grid [1,ntheta]
0043 %       - xpsin: normalized poloidal flux grid [1,nr]
0044 %       - psia_apRp: Edge normalized cylindrical poloidal flux (Wb) [1,1]
0045 %       - Xx: major radius [nr,ntheta]
0046 %       - Xy: elevation [nr,ntheta]
0047 %       - XBx: B-field component in R direction [nr,ntheta]
0048 %       - XBy: B-field component in Z direction [nr,ntheta]
0049 %       - XBPHI: B-field component in PHI direction [nr,ntheta]
0050 %       - xB0: minimum magnetic field on a flux-surface [1,nr]
0051 %       - xR0: corresponding major radius [1,nr]
0052 %       - xBT0: corresponding toroidal field [1,nr]
0053 %       - xBp0: corresponding poloidal field [1,nr]
0054 %       - xrho: normalized radius on the horizontal LFS midplane [1,nr]
0055 %       - xTe: electron temperature [1,nr]
0056 %       - xne: electron density [1,nr]
0057 %       - xzTi: ion temperatures [ni,nr]
0058 %       - xzni: ion densities [ni,nr]
0059 %       - zZi: ion charges [1,ni]
0060 %       - zmi: ion masses [1,ni]
0061 %       - xZeff: effective charge [1,nr]
0062 %       - Rp: major radius on axis [1,1]
0063 %       - Zp: elevation on axis [1,1]
0064 %       - Bax: magnetic field on axis [1,1]
0065 %       - ap: minor radius on LFS midplane [1,1]
0066 %       - Ip: Toroidal plasma current [1,1]
0067 %       - Ncoils: Number of toroidal magnetic field coils [1,1]
0068 %       - BeLi: Beta + li/2 [1,1]
0069 %       - xmhu0T2: pitch-angle (squared) at the trapped/passing boundary [1,nr]
0070 %       - xdV_dke: incremental FS volume at dke positions
0071 %       - xdA_dke: incremental FS poloidal surface at dke positions
0072 %       - HXR: Fast electron bremsstrahlung
0073 %
0074 %   By J. Decker RLE-MIT <jodecker@alum.mit.edu> and Y. Peysson DRFC-CEA <yves.peysson@cea.fr>
0075 %
0076 if nargin < 7,
0077     font = struct;
0078 end
0079 if nargin < 6,% individual figures are created
0080     flag_id = true;
0081     flag_consistency = true;
0082 else% graphs added to current object
0083     flag_consistency = false;
0084     if ~isobject(axs)
0085         flag_id = true;
0086         clear axs
0087     else
0088         flag_id = false;
0089     end
0090 end
0091 if nargin < 5
0092     mfactor = 10;%Factor for super theta grid
0093 end
0094 if nargin < 4
0095     method = 'linear';%Interpolation method
0096 end
0097 if nargin < 3
0098     display_mode = 0;
0099 elseif ~isstruct(display)
0100     display_mode = display;
0101 else
0102     display_mode = display.display_mode;
0103 end
0104 if nargin < 2,
0105     radialDKE = NaN;
0106 end
0107 if nargin < 1,
0108     error('Not enough input arguments in equilibrium_jd.m');
0109 end
0110 %
0111 if (~ischar(method) && isnan(method)) || isempty(method),
0112     method = 'linear';
0113 end
0114 %
0115 if isnan(mfactor) || isempty(mfactor),
0116     mfactor = 10;
0117 end
0118 %
0119 if isstruct(radialDKE),
0120     xpsin_S_dke = radialDKE.xpsin_S_dke;
0121     xpsin = radialDKE.xpsin_f;
0122 elseif length(radialDKE) > 1,
0123     xpsin_S_dke = radialDKE;
0124     xpsin = (xpsin_S_dke(1:end-1) + xpsin_S_dke(2:end))/2;
0125 else
0126     xpsin_S_dke = NaN;%no incremental volume or surface calculations
0127     xpsin = NaN;%the radial grid is the equilibrium grid
0128 end
0129 %
0130 % CASE OF LOCAL CALCULATION INDEPENDENT OF ANY ACTUAL TOKAMAK GENOMETRY
0131 %
0132 [npsi,ntheta] = size(equil.ptx);
0133 %
0134 Rp = equil.Rp;%Major radius of magnetic Axis
0135 Zp = equil.Zp;%Elevation of magnetic Axis
0136 Bax = equil.ptBPHI(1,1);%magnetic field on axis
0137 %
0138 if strcmp(equil.id,'LOCAL'),
0139     %
0140     % Note that circular flux-surfaces are assumed in this case
0141     %
0142     if npsi ~= 1 || ntheta ~= 1,
0143         error('LOCAL is reserved for equilibria with one (psi,theta) position only')
0144     end
0145     if isnan(xpsin)
0146         %disp('WARNING : xpsin is enforced to the value 0.5')
0147         xpsin = 0.5;
0148     end
0149     %
0150     theta = equil.theta;
0151     %
0152     xBT = abs(equil.ptBPHI);%Toroidal B-field
0153     xBp = sqrt(equil.ptBx.^2 + equil.ptBy.^2);%Poloidal B-field
0154     xB = sqrt(xBT.^2 + xBp.^2);%Total B-field
0155     %
0156     %index 0 refers to the minimum of the magnetic field
0157     %
0158     xr = sqrt(equil.ptx.^2 + equil.pty.^2);
0159     ep = xr/equil.Rp;
0160     xB0 = xB*(1 + ep*cos(theta))/(1 + ep);
0161     xBT0 = xBT*(1 + ep*cos(theta))/(1 + ep);%Toroidal field at B=Bmin
0162     xBp0 = xBp*(1 + ep*cos(theta))/(1 + ep);%Poroidal field at B=Bmin
0163     %
0164     xBmax = xB*(1 + ep*cos(theta))/(1 - ep);
0165     xmhu0T2 = 1 - xB0./xBmax; %Pitch-angle (squared) at the trapped-passing boundary
0166     %
0167     xx0 = xr;%
0168     %
0169     xrho = xpsin;
0170     xZeff = sum(equil.pzni.*equil.zZi.'.^2)./equil.pne;%Effective charge (a.u.)
0171     %
0172     ap = xr/xrho;
0173     xdV_2piRp_dke = NaN;
0174     xdA_dke = NaN;
0175     %
0176     equilDKE.tokname = equil.id;
0177     equilDKE.mtheta = theta;
0178     equilDKE.xpsin = xpsin;
0179     equilDKE.psia_apRp = equil.psi_apRp(end);
0180     equilDKE.Xx = equil.ptx;
0181     equilDKE.Xy = equil.pty;
0182     equilDKE.XBx = equil.ptBx;
0183     equilDKE.XBy = equil.ptBy;
0184     equilDKE.XBphi = equil.ptBPHI;
0185     equilDKE.xB0 = xB0;
0186     equilDKE.xx0 = xx0;
0187     equilDKE.xBT0 = xBT0;
0188     equilDKE.xBp0 = xBp0;
0189     equilDKE.xrho = xrho;
0190     equilDKE.xTe = equil.pTe;
0191     equilDKE.xne = equil.pne;
0192     equilDKE.xzTi = equil.pzTi;
0193     equilDKE.xzni = equil.pzni;
0194     equilDKE.zZi = equil.zZi;
0195     equilDKE.zmi = equil.zmi;
0196     equilDKE.xZeff = xZeff;
0197     equilDKE.Rp = equil.Rp;
0198     equilDKE.Zp = equil.Zp;
0199     equilDKE.Bax = Bax;
0200     equilDKE.ap = ap;
0201     equilDKE.xmhu0T2 = xmhu0T2;
0202     equilDKE.xdV_2piRp_dke = xdV_2piRp_dke;
0203     equilDKE.xdA_dke = xdA_dke;
0204     %
0205     return
0206 end
0207 %
0208 % FFT interpolation on R,Z, BR, BZ, BPHI as a function of theta (very precise, works for uniform theta grids only)
0209 %
0210 tpx = equil.ptx(:,1:(ntheta-1)).';
0211 tpy = equil.pty(:,1:(ntheta-1)).';
0212 tpBx = equil.ptBx(:,1:(ntheta-1)).';
0213 tpBy = equil.ptBy(:,1:(ntheta-1)).';
0214 tpBPHI = equil.ptBPHI(:,1:(ntheta-1)).';
0215 %
0216 if Rp < Inf,
0217     %
0218     % extrema removal
0219     %
0220     nc = (ntheta-1)/2 + 1;%central fft point
0221     %
0222     tol = 0.1;% tolerance on magnetic field variations
0223     %
0224     tpB0 = sqrt(tpBx.^2 + tpBy.^2 + tpBPHI.^2);% total magnetic field
0225     %
0226     opt_extr = 3;%1+1i;%
0227     %
0228     if real(opt_extr) == 1,
0229         %
0230         % Try cutting off high theta harmonics if there are secondary extrema
0231         %
0232         for t = 1:nc,
0233             %
0234             tpB = sqrt(tpBx.^2 + tpBy.^2 + tpBPHI.^2);% total magnetic field
0235             %
0236             var = max(max(abs(tpB - tpB0)./tpB0));
0237             %
0238             extr_ind = diff(sign(diff([tpB(end,:);tpB;tpB(1,:)]))) ~= 0;
0239             nextr = sum(extr_ind);
0240             %
0241             if t == nc || var > tol || all(nextr <= 2),
0242                 break
0243             end
0244             %
0245             ftpx = fft(tpx);
0246             ftpy = fft(tpy);
0247             ftpBx = fft(tpBx);
0248             ftpBy = fft(tpBy);
0249             ftpBPHI = fft(tpBPHI);
0250             %
0251             mask = (nc - t + 1):(nc + t - 1);
0252             %
0253             ftpx(mask,:) = 0;
0254             ftpy(mask,:) = 0;
0255             ftpBx(mask,:) = 0;
0256             ftpBy(mask,:) = 0;
0257             ftpBPHI(mask,:) = 0;
0258             %
0259             tpx = ifft(ftpx);
0260             tpy = ifft(ftpy);
0261             tpBx = ifft(ftpBx);
0262             tpBy = ifft(ftpBy);
0263             tpBPHI = ifft(ftpBPHI);
0264             %
0265         end
0266         %
0267         if t > 1,
0268             if var > tol || ~all(nextr <= 2),
0269                 %
0270                 disp(['Secondary extrema removal using Fourier filtering failed at  ',num2str(t - 1),'/',num2str(nc - 1),' harmonics.'])
0271                 %
0272                 if imag(opt_extr) == 1,
0273                     opt_extr = opt_extr + 1;
0274                 end
0275                 %
0276             else
0277                 disp(['Secondary extrema removal suceeded after cutting ',num2str(t - 1),'/',num2str(nc - 1),' harmonics.'])
0278             end
0279             %
0280             disp(['  -> variations on B : ',num2str(var),'/',num2str(tol),'.'])
0281         end
0282         %
0283     end
0284     %
0285     if real(opt_extr) == 2,        
0286         %
0287         % Try smoothing the poloidal field components to remove secondary extrema
0288         %
0289         tpx = equil.ptx(:,1:(ntheta-1)).';
0290         tpy = equil.pty(:,1:(ntheta-1)).';
0291         tpBx = equil.ptBx(:,1:(ntheta-1)).';
0292         tpBy = equil.ptBy(:,1:(ntheta-1)).';
0293         tpBPHI = equil.ptBPHI(:,1:(ntheta-1)).';
0294         %
0295         thetah = equil.theta(1:ntheta-1);
0296         dtheta = 2*pi/(ntheta-1);
0297         %
0298         for t = 1:nc,
0299             %
0300             tpB = sqrt(tpBx.^2 + tpBy.^2 + tpBPHI.^2);% total magnetic field
0301             %
0302             var = max(max(abs(tpB - tpB0)./tpB0));
0303             %
0304             extr_ind = diff(sign(diff([tpB(end,:);tpB;tpB(1,:)]))) ~= 0;
0305             nextr = sum(extr_ind);
0306             %
0307             if t == nc || var > tol || all(nextr <= 2),
0308                 break
0309             end
0310             %
0311             tpBx3 = smooth_jd([thetah-2*pi,thetah,thetah+2*pi],[tpBx;tpBx;tpBx],dtheta*t);
0312             tpBy3 = smooth_jd([thetah-2*pi,thetah,thetah+2*pi],[tpBy;tpBy;tpBy],dtheta*t);
0313             tpBx = tpBx3(ntheta:2*(ntheta-1),:);
0314             tpBy = tpBy3(ntheta:2*(ntheta-1),:);
0315             %
0316         end
0317         %
0318         if t > 1,
0319             if var > tol || ~all(nextr <= 2),
0320                 %
0321                 disp(['Secondary extrema removal using Gaussian smoothing failed at width ',num2str(dtheta*(t - 1)),'/pi.'])
0322                 %
0323                 if imag(opt_extr) == 1,
0324                     opt_extr = opt_extr + 1;
0325                 end
0326                 %
0327             else
0328                 disp(['Secondary extrema removal suceeded at width ',num2str(dtheta*(t - 1)),'/pi.'])
0329             end
0330             disp(['  -> variations on B : ',num2str(var),'/',num2str(tol),'.'])
0331         end
0332     end
0333     %
0334     if real(opt_extr) == 3,        
0335         %
0336         % Try removing secondary extrema directly, linear interpolation
0337         %
0338         tpx = equil.ptx(:,1:(ntheta-1)).';
0339         tpy = equil.pty(:,1:(ntheta-1)).';
0340         tpBx = equil.ptBx(:,1:(ntheta-1)).';
0341         tpBy = equil.ptBy(:,1:(ntheta-1)).';
0342         tpBPHI = equil.ptBPHI(:,1:(ntheta-1)).';
0343         %
0344         thetah = equil.theta(1:ntheta-1);
0345         %
0346         nt = 1000*nc;
0347         %
0348         for t = 1:nt,
0349             %
0350             tpB = sqrt(tpBx.^2 + tpBy.^2 + tpBPHI.^2);% total magnetic field
0351             %
0352             var = max(max(abs(tpB - tpB0)./tpB0));
0353             %
0354             extr_ind = diff(sign(diff([tpB(end,:);tpB;tpB(1,:)]))) ~= 0;%local extrema indices
0355             nextr = sum(extr_ind);
0356             %
0357             if t == nt || var > tol || all(nextr <= 2),
0358                 break
0359             end
0360             %
0361             disp(['Iteration : ',num2str(t),'; FS with extr. : ',num2str(sum(nextr > 2)),'/',...
0362                 num2str(npsi),'; B var. : ',num2str(var),'/',num2str(tol),'.'])
0363             %
0364             for ir = 1:npsi,
0365                 %
0366                 mask = extr_ind(:,ir);%extrema positions
0367                 %
0368                 itmax = find(tpB(:,ir) == max(tpB0(:,ir)),1,'first');
0369                 itmin = find(tpB(:,ir) == min(tpB0(:,ir)),1,'first');
0370                 %
0371                 mask(itmax) = false;%take global extrema indices away
0372                 mask(itmin) = false;%take global extrema indices away
0373                 %
0374                 if ~any(mask),
0375                     continue
0376                 end
0377                 %
0378                 tB = interp1([thetah(end) - 2*pi,thetah(~mask),2*pi],[tpB(end,ir);tpB(~mask,ir);tpB(1,ir)],thetah,'linear').';
0379                 %
0380                 tpBx(:,ir) = tpBx(:,ir).*tB./tpB(:,ir);
0381                 tpBy(:,ir) = tpBy(:,ir).*tB./tpB(:,ir);
0382                 tpBPHI(:,ir) = tpBPHI(:,ir).*tB./tpB(:,ir);
0383             end
0384             %
0385         end
0386         %
0387         if t > 1,
0388             if var > tol || ~all(nextr <= 2),
0389                 %
0390                 disp(['Secondary extrema direct removal failed after ',num2str(t-1),' iterations.'])
0391                 %
0392                 if imag(opt_extr) == 1,
0393                     opt_extr = opt_extr + 1;
0394                 end
0395                 %
0396             else
0397                 disp(['Secondary extrema direct removal suceeded after ',num2str(t-1),' iterations.'])
0398             end
0399             disp(['  -> variations on B : ',num2str(var),'/',num2str(tol),'.'])
0400         end
0401     end
0402     %
0403     if var > tol || ~all(nextr <= 2),
0404         error('Extrema removal failed, aborted')
0405     end    
0406     %
0407 end
0408 %
0409 for m = mfactor:-1:1,
0410     %
0411     mx = interpft(tpx,(ntheta-1)*m).'; %FFT interpolation
0412     my = interpft(tpy,(ntheta-1)*m).'; %FFT interpolation
0413     mBx = interpft(tpBx,(ntheta-1)*m).'; %FFT interpolation
0414     mBy = interpft(tpBy,(ntheta-1)*m).'; %FFT interpolation
0415     mBPHI = interpft(tpBPHI,(ntheta-1)*m).'; %FFT interpolation
0416     %
0417     if Rp < Inf
0418         %
0419         mB = sqrt(mBx.^2 + mBy.^2 + mBPHI.^2);
0420         %
0421         nextr = sum(diff(sign(diff(tpB))) ~= 0);
0422         %
0423         if all(nextr <= 2),
0424             break
0425         end
0426         %
0427     else
0428         break
0429     end
0430 end
0431 %
0432 mx = [mx,mx(:,1)]; %2*pi point retrieved
0433 my = [my,my(:,1)]; %2*pi point retrieved
0434 mBx = [mBx,mBx(:,1)]; %2*pi point retrieved
0435 mBy = [mBy,mBy(:,1)]; %2*pi point retrieved
0436 mBPHI = [mBPHI,mBPHI(:,1)]; %2*pi point retrieved
0437 %
0438 mtheta = 0:2*pi/((ntheta-1)*mfactor):2*pi; %super theta grid
0439 %
0440 if m < mfactor,
0441     %
0442     disp(['FT interpolation successful up to m = ',num2str(m),'/',num2str(mfactor)])
0443     %
0444     mmtheta = 0:2*pi/((ntheta-1)*m):2*pi; %super theta grid
0445     %
0446     mx = interp1(mmtheta,mx,mtheta);
0447     my = interp1(mmtheta,my,mtheta);
0448     mBx = interp1(mmtheta,mBx,mtheta);
0449     mBy = interp1(mmtheta,mBy,mtheta);
0450     mBPHI = interp1(mmtheta,mBPHI,mtheta);
0451     %
0452 end 
0453 %
0454 ap =  mx(npsi,1);%Plasma minor radius on the LFS midplane
0455 %
0456 psi_apRp = reshape(equil.psi_apRp,[1,npsi]);
0457 psia_apRp = psi_apRp(npsi);%Poloidal flux at the edge
0458 psin = psi_apRp/psia_apRp;%Normalized poloidal flux
0459 %
0460 if isfield(equil,'pzni') && isfield(equil,'pzTi') && isfield(equil,'zZi') && isfield(equil,'pne') && isfield(equil,'pTe') 
0461     pZeff = sum(equil.pzni.*repmat((equil.zZi.^2)',[1,npsi]))./equil.pne;%Effective charge (a.u.)
0462 end
0463 %
0464 if isnan(xpsin)
0465     xpsin = psin;
0466 end
0467 %
0468 nr = length(xpsin);
0469 %
0470 % Interpolation of calculation grid
0471 %
0472 Xx = interp1(psin,mx,xpsin,method);
0473 Xy = interp1(psin,my,xpsin,method);
0474 XBx = interp1(psin,mBx,xpsin,method);
0475 XBy = interp1(psin,mBy,xpsin,method);
0476 XBphi = interp1(psin,mBPHI,xpsin,method);
0477 %
0478 XBx(isnan(XBx)) = Inf;
0479 XBy(isnan(XBy)) = Inf;
0480 XBphi(isnan(XBphi)) = Inf;
0481 %
0482 xrho = Xx(:,1)'/ap;%Normalized radius on the horizontal LFS midplane
0483 %
0484 if isfield(equil,'pzni') && isfield(equil,'pzTi') && isfield(equil,'zZi') && isfield(equil,'pne') && isfield(equil,'pTe') 
0485     xTe = interp1(psin,equil.pTe,xpsin,method);
0486     xne = interp1(psin,equil.pne,xpsin,method);
0487     xzTi = interp1(psin,equil.pzTi.',xpsin,method).';
0488     xzni = interp1(psin,equil.pzni.',xpsin,method).';
0489     xZeff = interp1(psin,pZeff,xpsin,method);
0490 end
0491 %
0492 XBT = abs(XBphi);%Toroidal B-field
0493 XBp = sqrt(XBx.^2 + XBy.^2);%Poloidal B-field
0494 XB = sqrt(XBphi.^2 + XBx.^2 + XBy.^2);%Total B-field
0495 %
0496 % index 0 refers to the minimum of the magnetic field
0497 %
0498 [xB0,xiB0] = min(XB.'); %Minimum of B on flux surface
0499 xBmax = max(XB.'); %Maximum of B on flux surface
0500 %
0501 xmhu0T2 = 1 - xB0./xBmax; %Pitch-angle (squared) at the trapped-passing boundary
0502 %
0503 xiB0 = (xiB0 - 1)*nr + (1:nr);%Indices in (psi,theta) grid corresponding of B=Bmin
0504 xBT0 = XBT(xiB0);%Toroidal field at B=Bmin
0505 xBp0 = XBp(xiB0);%Poroidal field at B=Bmin
0506 xx0 = Xx(xiB0);%Major radius at B=Bmin
0507 %
0508 % volume and surface elements
0509 %
0510 % half grids in psi (for trapezoidal integration avoiding singularities)
0511 %
0512 hmx = (mx(1:npsi-1,:) + mx(2:npsi,:))/2;
0513 hmy = (my(1:npsi-1,:) + my(2:npsi,:))/2;
0514 hmBx = (mBx(1:npsi-1,:) + mBx(2:npsi,:))/2;
0515 hmBy = (mBy(1:npsi-1,:) + mBy(2:npsi,:))/2;
0516 hmBPHI = (mBPHI(1:npsi-1,:) + mBPHI(2:npsi,:))/2;
0517 %
0518 [hpq,hpqtilde,hpqbar,hpqhat,hpqcronos,hpB0] = qfactors_dke_yp(1,Rp,ap,hmx,hmy,hmBx,hmBy,hmBPHI);
0519 %
0520 dpsi_apRp = abs(diff(psi_apRp));
0521 %
0522 % FS Volume per 2pi*Rp
0523 %
0524 hpdV_2piRp = 2*pi*hpqhat./hpB0;
0525 pV_2piRp = [0,cumsum(hpdV_2piRp.*dpsi_apRp)];
0526 %
0527 % FS poloidal surface
0528 %
0529 hpdA = 2*pi*hpqbar./hpB0;
0530 pA = [0,cumsum(hpdA.*dpsi_apRp)];
0531 %
0532 % Other definitions for rho
0533 %
0534 psiT = 2*pi*[0,cumsum(hpq.*dpsi_apRp)];
0535 xrhoT = interp1(psin,sqrt(psiT/psiT(end)),xpsin);
0536 xrhoV = interp1(psin,sqrt(pV_2piRp/pV_2piRp(end)),xpsin);
0537 xrhoP = sqrt(xpsin);
0538 %
0539 % FS volume and poloidal surface
0540 %
0541 if ~isnan(xpsin_S_dke)
0542 %
0543 % DKE FS volume and poloidal surface
0544 %
0545     xV_2piRp_S_dke = interp1(psin,pV_2piRp,xpsin_S_dke,method);
0546     xA_S_dke = interp1(psin,pA,xpsin_S_dke,method);
0547 %
0548 % DKE FS invremental volume and poloidal surface
0549 %
0550     xdV_2piRp_dke = diff(xV_2piRp_S_dke);
0551     xdA_dke = diff(xA_S_dke);
0552 else
0553     xdV_2piRp_dke = NaN;
0554     xdA_dke = NaN;
0555 end
0556 %
0557 % Display of equilibrium properties
0558 %
0559 if real(display_mode) >= 2,
0560     %
0561     if real(display_mode) == 2,
0562         p_opt = input_dke_yp('Option for output figures : -1 (nothing), 0 (print), 1 (print and save), 2 (save)',-1,-1:2,'',[1,1]);
0563     else
0564         p_opt = -1;
0565     end
0566     %
0567     prho = mx(:,1)'/ap;
0568     %
0569     % displays the magnetic flux-surfaces geometry
0570     %
0571     if imag(display_mode) == 0,
0572         srho = linspace(0,1,21);
0573     else
0574         srho = linspace(0,1,11);
0575     end
0576     %
0577     nthetad = 25;
0578     stheta = linspace(0,2*pi,nthetad);
0579     sx = interp1(prho,mx,srho,method);  
0580     sy = interp1(prho,my,srho,method); 
0581     stx = interp1(mtheta,mx',stheta,method)';  
0582     sty = interp1(mtheta,my',stheta,method)';    
0583     %
0584     xmin = min(min(equil.ptx));
0585     xmax = max(max(equil.ptx));
0586     ymin = min(min(equil.pty));
0587     ymax = max(max(equil.pty));
0588     %
0589     if flag_id,
0590         figure('Name','Magnetic equilibrium flux surface contour'),clf
0591         ax = gca;
0592         red = 0.9;
0593         lspace = 0.7;
0594         bspace = 0.5;
0595         tit = 'Flux-surfaces labeled by \rho=0,0.05,0.1...1';
0596     else
0597         ax = axs(1);
0598         red = 1;
0599         lspace = NaN;
0600         bspace = NaN;
0601         tit = '';
0602     end
0603     %
0604     if imag(display_mode) == 0,
0605         [ax] = graph1D_jd(sx',sy',0,0,'x(m)','y(m)','',NaN,[xmin,xmax],[ymin,ymax],'-','none','r',2,font,ax,red,lspace,bspace);drawnow
0606         [ax] = graph1D_jd(stx,sty,0,0,'','',tit,NaN,[xmin,xmax],[ymin,ymax],'-','none','b',0.5,font,ax);
0607     else
0608         [ax] = graph1D_jd(sx',sy',0,0,'x(m)','y(m)','',NaN,[xmin,xmax],[ymin,ymax],'-','none','r',0.5,font,ax,red,lspace,bspace);drawnow
0609     end
0610     axis(ax,'equal');drawnow  
0611     print_jd(p_opt,['Fig_EQUIL_',equil.id,'_mag'])
0612     %
0613     if isfield(equil,'pzni') && isfield(equil,'pzTi') && isfield(equil,'zZi') && isfield(equil,'pne') && isfield(equil,'pTe') 
0614         if flag_id,
0615             figure('Name','Electron Temperature and Density Profiles'),clf
0616             titlefig = '';
0617             if isfield(equil,'Te_vol'), titlefig = ['<T_e>=',num2str(round(equil.Te_vol*100)/100),'keV'];end
0618             if isfield(equil,'ne_vol'), titlefig = [titlefig,', <n_e>=',num2str(round(100*equil.ne_vol/1e19)/100),'x10^{+19} m^{-3}'];end
0619             [ax] = graph1D2jd(prho,xTe,prho,xne/1e19,'\rho','T_e (keV)','n_e (10^{19} m^{-3})',titlefig,...
0620                 [0,1],NaN,NaN,'-','none','r',2,'--','none','b',2,20,0.9,0.5,0.7);
0621             set(ax(1),'xtick',[0:0.2:1]);set(ax(2),'xtick',[0:0.2:1])
0622             %
0623             print_jd(p_opt,['Fig_EQUIL_',equil.id,'_tene'])
0624         elseif length(axs) >= 2,
0625             if imag(display_mode) == 0,
0626                 ax = axs(2);
0627                 leg = {'T_e (keV)','n_e (10^{19} m^{-3})','Location','NorthEast'};
0628                 [ax] = graph1D_jd(prho,equil.pTe,0,0,'','','',NaN,[0,1],NaN,'-','none','r',2,font,ax);
0629                 [ax] = graph1D_jd(prho,equil.pne/1e19,0,0,'\rho','','',leg,[0,1],NaN,'--','none','b',2,font,ax);
0630                 set(ax,'xtick',[0:0.2:1]);
0631             else
0632                 ax = axs(2);
0633                 %
0634                 nphi = 361;
0635                 nm = length(mtheta);
0636                 sR = Rp + [sx(end:-1:2,(nm + 1)/2);sx(:,1)];
0637                 sphi = linspace(0,2*pi,nphi);
0638                 XR = repmat(sR.',[nphi,1]);
0639                 Xphi = repmat(sphi.',[1,21]);
0640                 Rmax = max(sR);
0641                 %
0642                 [ax] = graph1D_jd(XR.*cos(Xphi),XR.*sin(Xphi),0,0,'X(m)','Y(m)','',NaN,[-Rmax,Rmax],[-Rmax,Rmax],'-','none','r',0.5,font,ax,red,lspace,bspace);
0643                 axis(ax,'equal');set(ax,'xtickmode','auto','ytickmode','auto');drawnow;                
0644             end
0645         end
0646     end
0647     %
0648     % displays the equilibrium profiles
0649     %
0650     if isfield(equil,'pzni') && isfield(equil,'pzTi') && isfield(equil,'zZi') && isfield(equil,'pne') && isfield(equil,'pTe') 
0651         %
0652         if flag_id, 
0653             figure('Name','Equilibrium profile parameters'),clf
0654             set(gcf,'units','inches','position',[0.5,1,10,6]);
0655             set(gca,'position',[0,0,1,1]);
0656             %
0657             subplot('position',[1/12,0.5+1/16,5/24,5/16])
0658               [ax] = graph1D_jd(prho,psin,0,0,'','','',NaN,[0,1],[-eps,1],'-','none','b',0.5);
0659               [ax] = graph1D_jd(xrho,xpsin,0,0,'\rho','\psi_n','',NaN,[0,1],[-eps,1],'none','+','r',2);
0660             subplot('position',[1/3+1/12,0.5+1/16,5/24,5/16])
0661               [ax] = graph1D_jd(prho,equil.pTe,0,0,'','','',NaN,[0,1],[-eps,max(equil.pTe)],'-','none','b',0.5);
0662               [ax] = graph1D_jd(xrho,xTe,0,0,'\rho','T_e (keV)',equil.id,NaN,[0,1],[-eps,max(equil.pTe)],'none','+','r',2);
0663             subplot('position',[2/3+1/12,0.5+1/16,5/24,5/16])
0664               [ax] = graph1D_jd(prho,equil.pne/1e19,0,0,'','','',NaN,[0,1],[-eps,max(equil.pne/1e19)],'-','none','b',0.5);
0665               [ax] = graph1D_jd(xrho,xne/1e19,0,0,'\rho','n_e (10^{19} m^{-3})','',NaN,[0,1],[-eps,max(equil.pne/1e19)],'none','+','r',2);
0666             subplot('position',[1/12,1/8,5/24,5/16])
0667               [ax] = graph1D_jd(prho,equil.pzTi,0,0,'','','',NaN,[0,1],[-eps,max(max(equil.pzTi))],'-','none','b',0.5);
0668               [ax] = graph1D_jd(xrho,xzTi,0,0,'\rho','T_i (keV)','',NaN,[0,1],[-eps,max(max(equil.pzTi))],'none','+','r',2);
0669             subplot('position',[1/3+1/12,1/8,5/24,5/16])
0670               [ax] = graph1D_jd(prho,equil.pzni/1e19,0,0,'','','',NaN,[0,1],[-eps,max(max(equil.pzni/1e19))],'-','none','b',0.5);
0671               [ax] = graph1D_jd(xrho,xzni/1e19,0,0,'\rho','n_i (10^{19} m^{-3})','',NaN,[0,1],[-eps,max(max(equil.pzni/1e19))],'none','+','r',2);
0672             subplot('position',[2/3+1/12,1/8,5/24,5/16])
0673               [ax] = graph1D_jd(prho,pZeff,0,0,'','','',NaN,[0,1],[-eps,max(pZeff)],'-','none','b',0.5);
0674               [ax] = graph1D_jd(xrho,xZeff,0,0,'\rho','Z_{eff}','',NaN,[0,1],[-eps,max(pZeff)],'none','+','r',2);
0675         elseif length(axs) >= 3,
0676             ax = axs(3);
0677             [ax] = graph1D_jd(prho,pZeff,0,0,'\rho','Z_{eff}','',NaN,[0,1],[0,ceil(max(pZeff))],'-','none','r',2,font,ax,red,lspace,bspace);drawnow     
0678             set(ax,'ytick',[0:ceil(max(pZeff))]);
0679         end
0680         %
0681         if p_opt >= 0,
0682             %
0683             axlist = get(gcf,'children');
0684             figure('Name','Equilibrium profile > print'),clf
0685             cpaxlist = copyobj(axlist,gcf);
0686             %
0687             set(gcf,'position',[0,0,600,900]),
0688             set(cpaxlist(1),'units','normalized','position',[0.65,0.1,0.3,0.2])
0689             set(cpaxlist(2),'units','normalized','position',[0.65,0.4,0.3,0.2])
0690             set(cpaxlist(3),'units','normalized','position',[0.15,0.4,0.3,0.2])
0691             set(cpaxlist(4),'units','normalized','position',[0.65,0.7,0.3,0.2])
0692             set(cpaxlist(5),'units','normalized','position',[0.15,0.7,0.3,0.2])
0693             set(cpaxlist(6),'units','normalized','position',[0.15,0.1,0.3,0.2])    
0694             %
0695             set(get(cpaxlist(5),'title'),'string','')
0696             %
0697             print_jd(p_opt,['Fig_EQUIL_',equil.id,'_prof']);    
0698             %
0699         end
0700     end
0701 end
0702 %
0703 if isfield(equil,'pzni') && isfield(equil,'pzTi') && isfield(equil,'zZi') && isfield(equil,'pne') && isfield(equil,'pTe') 
0704     zZi = equil.zZi;
0705     zmi = equil.zmi;
0706 end
0707 %
0708 if real(display_mode) >=1,
0709     disp('-------------------------------------------------------------------------------------------------------------------');
0710     disp(['- Scenario name = ',equil.id]);
0711     disp('-------------------------------------------------------------------------------------------------------------------');
0712     disp(['------------------------------------- Tokamak magnetic equilibrium parameters ------------------------------------']);
0713     disp(['- Plasma minor radius on LFS midplane (m) ap = ',num2str(ap)]);
0714     disp(['- Plasma major radius on axis (m) Rp = ',num2str(Rp)]);
0715     disp(['- Plasma elevation on axis (m) Zp = ',num2str(Zp)]);
0716     disp(['- Edge normalized cylindrical poloidal flux (Wb) = ',num2str(psia_apRp)]);
0717     if isfield(equil,'signs') && ~isempty(equil.signs),
0718         disp(['- Sign of the plasma current (LUKE convention) = ',num2str(equil.signs.Ip)]);
0719         disp(['- Sign of the toroidal magnetic field (LUKE convention) = ',num2str(equil.signs.Bt)]);
0720     end
0721     disp('-------------------------------------------------------------------------------------------------------------------');
0722 end
0723 %
0724 if flag_consistency && real(display_mode) >= 1,
0725     disp('Calculation of the equilibrium consistency (wait few seconds).');
0726     equil = equilconsistency_yp(equil,'',real(display_mode));
0727 end
0728 %
0729 equilDKE.tokname = equil.id;
0730 equilDKE.mtheta = mtheta;
0731 equilDKE.xpsin = xpsin;
0732 equilDKE.psia_apRp = psia_apRp;
0733 equilDKE.Xx = Xx;
0734 equilDKE.Xy = Xy;
0735 equilDKE.XBx = XBx;
0736 equilDKE.XBy = XBy;
0737 equilDKE.XBphi = XBphi;
0738 equilDKE.xB0 = xB0;
0739 equilDKE.xx0 = xx0;
0740 equilDKE.xBT0 = xBT0;
0741 equilDKE.xBp0 = xBp0;
0742 equilDKE.xrho = xrho;
0743 equilDKE.xrhoP = xrhoP;
0744 equilDKE.xrhoT = xrhoT;
0745 equilDKE.xrhoV = xrhoV;
0746 %
0747 if isfield(equil,'pzni') && isfield(equil,'pzTi') && isfield(equil,'zZi') && isfield(equil,'pne') && isfield(equil,'pTe') 
0748     equilDKE.xTe = xTe;
0749     equilDKE.xne = xne;
0750     equilDKE.xzTi = xzTi;
0751     equilDKE.xzni = xzni;
0752     equilDKE.zZi = zZi;
0753     equilDKE.zmi = zmi;
0754     equilDKE.xZeff = xZeff;
0755 end
0756 %
0757 equilDKE.Rp = Rp;
0758 equilDKE.Zp = Zp;
0759 equilDKE.Bax = Bax;
0760 equilDKE.ap = ap;
0761 equilDKE.xmhu0T2 = xmhu0T2;
0762 equilDKE.xdV_2piRp_dke = xdV_2piRp_dke;
0763 equilDKE.xdA_dke = xdA_dke;
0764 %
0765 if isfield(equil,'equil_fit'),
0766     equilDKE.equil_fit = equil.equil_fit;
0767 end
0768 %
0769 if isfield(equil,'ip'),
0770     equilDKE.ip = equil.ip;
0771 end
0772 %
0773 if isfield(equil,'li'),
0774     equilDKE.li = equil.li;
0775 end
0776 %
0777 if nargout >= 1,
0778     varargout{1} = equilDKE;
0779 end
0780 
0781
Community support and wiki are available on Redmine. Last update: 18-Apr-2019.