0001 function [equil,gradpthermalsgradrho,gradpmagsgradrho,gradpkinsgradrho,jparz_fsav,jparz_ohm_fsav] = equilconsistency_yp(equil,equil_fit,display_mode,prho,mtheta)
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024
0025
0026
0027
0028 if nargin == 1,
0029 if isfield(equil,'equil_fit'),
0030 equil_fit = equil.equil_fit;
0031 else
0032 info_dke_yp(2,['Vectorized form of the magnetic equilibrium ',equil.id,' is being calculated...'],0);
0033 equil_fit = fitequil_yp(equil);
0034 info_dke_yp(2,' ...done',0);
0035 end
0036
0037 display_mode = 0;
0038 prho = max(101,length(equil.psi_apRp));
0039 mtheta = max(105,length(equil.theta));
0040 end
0041
0042 if nargin == 2,
0043 display_mode = 0;
0044 prho = max(101,length(equil.psi_apRp));
0045 mtheta = max(105,length(equil.theta));
0046 end
0047
0048 if nargin == 3,
0049 prho = max(101,length(equil.psi_apRp));
0050 mtheta = max(105,length(equil.theta));
0051 end
0052
0053 if nargin == 4,
0054 mtheta = max(105,length(equil.theta));
0055 end
0056
0057 rho = linspace(0.01,1,prho);
0058 drho = rho(2)-rho(1);
0059
0060 theta = linspace(0,2*pi,mtheta);
0061 dtheta = theta(2)-theta(1);
0062
0063 if isempty(equil_fit),
0064 if isfield(equil,'equil_fit'),
0065 equil_fit = equil.equil_fit;
0066 else
0067 info_dke_yp(2,['Vectorized form of the magnetic equilibrium ',equil.id,' is being calculated...'],0);
0068 equil_fit = fitequil_yp(equil);
0069 info_dke_yp(2,' ...done',0);
0070 equil.equil_fit = equil_fit;
0071 end
0072 end
0073
0074 equilval = equilval_yp(equil_fit,rho,theta);
0075
0076 Te = equilval.Te;
0077 dTedrho = equilval.dTedrho;
0078 d2Tedrho2 = equilval.d2Tedrho2;
0079
0080 ne = equilval.ne;
0081 dnedrho = equilval.dnedrho;
0082 d2nedrho2 = equilval.d2nedrho2;
0083
0084 zTi = equilval.zTi;
0085 dzTidrho = equilval.dzTidrho;
0086 d2zTidrho2 = equilval.d2zTidrho2;
0087
0088 zni = equilval.zni;
0089 dznidrho = equilval.dznidrho;
0090 d2znidrho2 = equilval.d2znidrho2;
0091
0092 x = equilval.x;
0093 dxdtheta = equilval.dxdtheta;
0094 d2xdtheta2 = equilval.d2xdtheta2;
0095 dxdrho = equilval.dxdrho;
0096 d2xdthetadrho = equilval.d2xdthetadrho;
0097
0098 y = equilval.y;
0099 dydtheta = equilval.dydtheta;
0100 d2ydtheta2 = equilval.d2ydtheta2;
0101 dydrho = equilval.dydrho;
0102 d2ydthetadrho = equilval.d2ydthetadrho;
0103
0104 r = equilval.r;
0105 drdtheta = equilval.drdtheta;
0106 drdrho = equilval.drdrho;
0107
0108 salpha = equilval.salpha;
0109 dsalphadtheta = equilval.dsalphadtheta;
0110 dsalphadrho = equilval.dsalphadrho;
0111
0112 calpha = equilval.calpha;
0113 dcalphadtheta = equilval.dcalphadtheta;
0114 dcalphadrho = equilval.dcalphadrho;
0115
0116 Bx = equilval.Bx;
0117 dBxdtheta = equilval.dBxdtheta;
0118 d2Bxdtheta2 = equilval.d2Bxdtheta2;
0119 dBxdrho = equilval.dBxdrho;
0120 d2Bxdthetadrho = equilval.d2Bxdthetadrho;
0121
0122 By = equilval.By;
0123 dBydtheta = equilval.dBydtheta;
0124 d2Bydtheta2 = equilval.d2Bydtheta2;
0125 dBydrho = equilval.dBydrho;
0126 d2Bydthetadrho = equilval.d2Bydthetadrho;
0127
0128 Bz = equilval.Bz;
0129 dBzdtheta = equilval.dBzdtheta;
0130 d2Bzdtheta2 = equilval.d2Bzdtheta2;
0131 dBzdrho = equilval.dBzdrho;
0132 d2Bzdthetadrho = equilval.d2Bzdthetadrho;
0133
0134 BP = equilval.BP;
0135 dBPdtheta = equilval.dBPdtheta;
0136 dBPdrho = equilval.dBPdrho;
0137
0138 B = equilval.B;
0139 dBdtheta = equilval.dBdtheta;
0140 dBdrho = equilval.dBdrho;
0141
0142 gradrho = equilval.gradrho;
0143 dgradrhodtheta = equilval.dgradrhodtheta;
0144 dgradrhodrho = equilval.dgradrhodrho;
0145
0146 Upsilon = equilval.Upsilon;
0147 dUpsilondrho = equilval.dUpsilondrho;
0148 dUpsilondtheta = equilval.dUpsilondtheta;
0149
0150 psin = equilval.psin;
0151 dpsindrho = equilval.dpsindrho;
0152 dxstardpsin = equilval.dxstardpsin;
0153 ddxstardpsindrho = equilval.ddxstardpsindrho;
0154
0155 [qe,me,mp,mn,e0,mu0,re,mc2,clum,alpha,kB] = pc_dke_yp;
0156
0157 if isfield(equil,'eqdsk')
0158 equil.eqdsk.prho = psi2rho_jd(equil,equil.eqdsk.psin,'spline');
0159 equil.eqdsk.gradpkinsgradrho = gradient(equil.eqdsk.pp)./gradient(equil.eqdsk.prho);
0160 end
0161
0162 jparz_fsav = zeros(prho,1);
0163 jparz_ohm_fsav = zeros(prho,1);
0164 jperpz_fsav = zeros(prho,1);
0165
0166 Ip = 0;
0167 Ip_ohm = 0;
0168
0169 S = zeros(prho,1);
0170 Surf = 0;
0171
0172 if ~isempty(zni) && ~isempty(zTi) && ~isempty(Te) && ~isempty(ne')
0173 pthermal_ion = sum(zTi.*zni,2);
0174 dpthermal_iondrho = sum(zTi.*dznidrho + dzTidrho.*zni,2);
0175
0176 pthermal = (Te.*ne + pthermal_ion)*qe*1000;
0177 gradpthermalsgradrho = (Te.*dnedrho + dTedrho.*ne + dpthermal_iondrho)*qe*1000;
0178 else
0179 pthermal_ion = NaN;
0180 dpthermal_iondrho = NaN;
0181 pthermal = NaN;
0182 gradpthermalsgradrho = NaN;
0183 end
0184
0185 a1 = drdrho.*BP./calpha + r.*dBPdrho./calpha - r.*BP.*dcalphadrho./calpha./calpha;
0186 a2 = dsalphadtheta./Upsilon./calpha - salpha.*dcalphadtheta./Upsilon./calpha./calpha - salpha.*dUpsilondtheta./calpha./Upsilon./Upsilon;
0187 a4 = dUpsilondrho.*Bz + Upsilon.*dBzdrho;
0188
0189 b1 = -a1.*BP.*calpha./r;
0190 b2 = a2.*BP.*BP.*Upsilon.*calpha./r./gradrho;
0191 b4 = -a4.*Bz./Upsilon;
0192
0193 c1 = a1.*Bz.*gradrho.*calpha./r./B;
0194 c2 = -a2.*BP.*Upsilon.*calpha.*Bz./r./B;
0195 c4 = -a4.*BP.*gradrho./Upsilon./B;
0196
0197 d1 = a1.*gradrho.*calpha./r;
0198 d2 = -a2.*BP.*Upsilon.*calpha./r;
0199
0200 e1 = a1.*BP.*BP.*gradrho.*calpha./r./B./B;
0201 e2 = -a2.*BP.*BP.*BP.*Upsilon.*calpha./r./B./B;
0202 e4 = a4.*BP.*Bz.*gradrho./Upsilon./B./B;
0203
0204 gradpmagsgradrho = (b1 + b2 + b4)./mu0;
0205
0206 jz = (d1 + d2)./mu0;
0207
0208 jpar = (c1 + c2 + c4)./mu0;
0209 jparz = jpar.*Bz./B;
0210
0211 jperpz = (e1 + e2 + e4)./mu0;
0212
0213 S = sum(r./Upsilon./BP./calpha.*dtheta,2);
0214 Surf = sum(sum(abs(r./gradrho./calpha)*dtheta*drho));
0215
0216 jz_fsav = sum(jz.*r./Upsilon./BP./calpha.*dtheta,2)./S;
0217 jparz_fsav = sum(jparz.*r./Upsilon./BP./calpha.*dtheta,2)./S;
0218 jperpz_fsav = sum(jperpz.*r./Upsilon./BP./calpha.*dtheta,2)./S;
0219
0220 Ip_z = sum(sum(jz.*r.*dtheta*drho./gradrho./calpha));
0221 Ip_z_fsav = abs(equil.psi_apRp(end))*sum(jz_fsav.*S./dxstardpsin(:)*drho);
0222
0223 Ip_parperpz = sum(sum((jparz + jperpz).*r.*dtheta*drho./gradrho./calpha));
0224 Ip_parperpz_fsav = abs(equil.psi_apRp(end))*sum((jparz_fsav + jperpz_fsav).*S./dxstardpsin(:)*drho);
0225
0226 F = Bz.*(x + equil.Rp);
0227 dFdpsi = (dBzdrho.*(x + equil.Rp) + Bz.*dxdrho).*(dxstardpsin(:)*ones(1,mtheta))/equil.psi_apRp(end)*equil.Rp;
0228 FdFdpsi = F.*dFdpsi;
0229
0230 if isfield(equil,'eqdsk')
0231
0232 equilDKE = equilibrium_jd(equil,equil.eqdsk.psin,0,'spline');
0233
0234 rhoP = interp1([0,equilDKE.xrho,1],[0,equilDKE.xrhoP,1],rho);
0235 rhoT = interp1([0,equilDKE.xrho,1],[0,equilDKE.xrhoT,1],rho);
0236 rhoV = interp1([0,equilDKE.xrho,1],[0,equilDKE.xrhoV,1],rho);
0237
0238 xq = qfactors_dke_yp(1,equilDKE.Rp,equilDKE.ap,equilDKE.Xx,equilDKE.Xy,equilDKE.XBx,equilDKE.XBy,equilDKE.XBphi,equilDKE.xx0,equilDKE.xB0,equilDKE.xBT0,equilDKE.xBp0);
0239
0240 q = equilDKE.ap/equilDKE.Rp*interp1(equilDKE.xrho,xq,rho);
0241
0242 pF_eqdsk = spline(equil.eqdsk.prho,equil.eqdsk.pF,rho)'*ones(1,length(theta));
0243 pFdFdpsi_eqdsk = spline(equil.eqdsk.prho,equil.eqdsk.pFdFdpsi,rho)'*ones(1,length(theta));
0244 figure,plot(rho,pFdFdpsi_eqdsk,'r',rho,FdFdpsi/3.5,'b')
0245
0246 pp_eqdsk = spline(equil.eqdsk.prho,equil.eqdsk.pp,rho)'*ones(1,length(theta));
0247 pdpdpsi_eqdsk = spline(equil.eqdsk.prho,equil.eqdsk.pdpdpsi,rho)'*ones(1,length(theta));
0248 pdpdpsi_eqdsk1 = gradient(equil.eqdsk.pp)./gradient(equil.eqdsk.psin)*abs(equil.eqdsk.psia);
0249
0250
0251 gradpkinsgradrho = equil.eqdsk.gradpkinsgradrho;
0252
0253
0254
0255
0256 jz_eqdsk = Upsilon.*pdpdpsi_eqdsk*equilDKE.Rp + pFdFdpsi_eqdsk./Upsilon/equilDKE.Rp;
0257
0258
0259 if isfield( equil.eqdsk,'Ip'),
0260 Ip_eqdsk = equil.eqdsk.Ip;
0261 end
0262 else
0263 gradpkinsgradrho = [];
0264 jphi_eqdsk = [];
0265 Ip_eqdsk = [];
0266 end
0267
0268 for j = 1:6,
0269 jparz_fsav_cyl(:,j) = j*abs(Ip_z)*(1 - rho(:).^2).^(j-1)/pi/equil_fit.ap/equil_fit.ap/1e6;
0270 end
0271
0272 Surf_cyl = pi*equil_fit.ap*equil_fit.ap;
0273
0274 if ~isempty(zni) && ~isempty(zTi) && ~isempty(Te) && ~isempty(ne),
0275 jparz_ohm = (Te*ones(1,mtheta)).^1.5;
0276 Ip_parz_ohm= sum(sum(jparz_ohm.*r.*dtheta*drho./gradrho./calpha));
0277
0278 jparz_ohm_fsav = sum(jparz_ohm.*r./Upsilon./BP./calpha.*dtheta,2)./S;
0279 Ip_parz_ohm_fsav = abs(equil.psi_apRp(end))*sum(jparz_ohm_fsav.*S./dxstardpsin(:)*drho);
0280
0281 jparz_ohm_fsav = abs(jparz_ohm_fsav*Ip_parperpz_fsav/Ip_parz_ohm_fsav./S).*sign(jparz_fsav);
0282 else
0283 jparz_ohm = NaN;
0284 jparz_ohm_fsav = NaN;
0285 Ip_parz_ohm_fsav = NaN;
0286 end
0287
0288 L = sum(r(end,:).*dtheta./calpha(end,:));
0289 LB = sum(r(end,:).*dtheta.*BP(end,:)./calpha(end,:));
0290
0291 Ip_contour = LB/mu0;
0292
0293 VB2 = sum(sum(r.*Upsilon.*drho.*dtheta.*BP.*BP./gradrho./calpha));
0294 V = sum(sum(r.*Upsilon.*drho.*dtheta./gradrho./calpha));
0295
0296 XTe = Te(:)*ones(1,mtheta);
0297 VTe = sum(sum(r.*Upsilon.*drho.*dtheta.*XTe./gradrho./calpha));
0298 Te_vol = VTe/V;
0299
0300 Xne = ne(:)*ones(1,mtheta);
0301 Vne = sum(sum(r.*Upsilon.*drho.*dtheta.*Xne./gradrho./calpha));
0302 ne_vol = Vne/V;
0303
0304 for iz = 1:length(equil.zZi)
0305 XTi = zTi(:,iz)*ones(1,mtheta);
0306 VTi = sum(sum(r.*Upsilon.*drho.*dtheta.*XTi./gradrho./calpha));
0307 zTi_vol(iz) = VTi/V;
0308
0309 Xni = zni(:,iz)*ones(1,mtheta);
0310 Vni = sum(sum(r.*Upsilon.*drho.*dtheta.*Xni./gradrho./calpha));
0311 zni_vol(iz) = Vni/V;
0312 end
0313
0314 li = L*L*VB2/LB/LB/V;
0315
0316 [pshift,pelong,ptrian,pli] = equilmoments_jd(equil.theta,equil.ptx,equil.pty,mtheta);
0317
0318 disp(['Plasma current, Ip(from j_z) = ',num2str(Ip_z/1e6),' (MA)']);
0319 disp(['Plasma current, Ip(from j_z_fsav) = ',num2str(Ip_z_fsav/1e6),' (MA)']);
0320 disp(['Plasma current, Ip(from j_par_z + j_perp_z) = ',num2str(Ip_parperpz/1e6),' (MA)']);
0321 disp(['Plasma current, Ip(from j_par_z_fsav + j_perp_z_fsav) = ',num2str(Ip_parperpz_fsav/1e6),' (MA)']);
0322 disp(['Plasma current, Ip(from contour) = ',num2str(Ip_contour/1e6),' (MA)']);
0323 if ~isempty(Ip_eqdsk),
0324 disp(['Plasma current, Ip(eqdsk) = ',num2str(Ip_eqdsk/1e6),' (MA)']);
0325 end
0326 disp('-------------------------------------------------------------')
0327 disp(['Normalized internal inductance, li = ',num2str(li)]);
0328 disp('-------------------------------------------------------------')
0329 disp(['Shafranov shift = ',num2str(pshift(end))]);
0330 disp(['Plasma elongation = ',num2str(pelong(end))]);
0331 disp(['Plasma triangularity = ',num2str(ptrian(end))]);
0332
0333 if equil.psi_apRp(end)*equil.ptBy(end,1) < 0,
0334 disp('******** WARNING : equilibrium signs not self-consistent ********')
0335 if equil.psi_apRp(end) > 0,
0336 disp('******** Poloidal magnetic flux corresponds to positive plasma current (clockwise when tokamak viewed from the top) ********')
0337 disp('******** Poloidal magnetic field direction is clockwise corresponding to a negative plasma current (counter-clockwise when tokamak viewed from the top) ********')
0338 else
0339 disp('******** Poloidal magnetic flux corresponds to negative plasma current (counter-clockwise when tokamak viewed from the top) ********')
0340 disp('******** Poloidal magnetic field direction is counter-clockwise corresponding to a positive plasma current (clockwise when tokamak viewed from the top) ********')
0341 end
0342
0343 disp('******** The problem must be fixed before any C3PO/LUKE simulation ********')
0344 disp('')
0345 opt_gui = input_dke_yp('Do you want to fix the problem ? (yes:1,no:0)',0,[0,1]) == 1;
0346 if opt_gui == 1,
0347 opt_gui1 = input_dke_yp('Is the plasma current clockwise (1) or anticlockwise (-1) when tokamak viewed from the top ?',1,[-1,1]) == 1;
0348 if opt_gui1 == 1,
0349 if equil.psi_apRp(end) > 0
0350 equil.ptBx = - equil.ptBx;
0351 equil.ptBy = - equil.ptBy;
0352 else
0353 equil.psi_apRp = -1*equil.psi_apRp;
0354 end
0355 elseif opt_gui1 == -1,
0356 if equil.psi_apRp(end) > 0
0357 equil.psi_apRp = -1*equil.psi_apRp;
0358 else
0359 equil.ptBx = - equil.ptBx;
0360 equil.ptBy = - equil.ptBy;
0361 end
0362 end
0363
0364 end
0365 end
0366
0367 if equil.psi_apRp(end) > 0,
0368 disp('-------> Poloidal magnetic field is self-consistent with Ip positive (clockwise when tokamak viewed from the top)')
0369 elseif equil.psi_apRp(end) < 0
0370 disp('-------> Poloidal magnetic field is self-consistent with Ip negative (counter-clockwise when tokamak viewed from the top)')
0371 end
0372
0373 if equil.ptBPHI(1,1) > 0
0374 disp('-------> Toroidal magnetic field is positive (clockwise when tokamak viewed from the top)')
0375 else
0376 disp('-------> Toroidal magnetic field is negative (counter-clockwise when tokamak viewed from the top)')
0377 end
0378
0379 if display_mode == 2,
0380
0381 figure('Name','Pressure gradient'),clf
0382 if isfield(equil,'eqdsk')
0383 if ~isnan(jparz_ohm_fsav),
0384 leg = {'Thermal','Kinetic','Magnetic','Location','SouthWest'};
0385 else
0386 leg = {'Kinetic','Magnetic','Location','SouthWest'};
0387 end
0388 else
0389 if ~isnan(jparz_ohm_fsav),
0390 leg = {'Thermal','Magnetic','Location','SouthWest'};
0391 else
0392 leg = {'Magnetic','Location','SouthWest'};
0393 end
0394 end
0395 if ~isnan(jparz_ohm_fsav),
0396 [ax] = graph1D_jd(rho,gradpthermalsgradrho(:,1),0,0,'\rho','\nablap\cdot\psi/|\nabla\rho|',strrep(equil.id,'_','-'),NaN,NaN,NaN,'-','none','r',2);drawnow
0397 end
0398 if isfield(equil,'eqdsk')
0399 [ax] = graph1D_jd(equil.eqdsk.prho,equil.eqdsk.gradpkinsgradrho,0,0,'\rho','\nablap\cdot\psi/|\nabla\rho|',strrep(equil.id,'_','-'),NaN,NaN,NaN,'-','none','g',1);
0400 [ax] = graph1D_jd(rho,gradpmagsgradrho,0,0,'\rho','\nablap\cdot\psi/|\nabla\rho|',[strrep(equil.id,'_','-'),', \theta\in[0:\pi/',int2str((mtheta-1)/2),':2\pi]'],NaN,NaN,NaN,'--','none','b',2);
0401 if ~isnan(jparz_ohm_fsav),
0402 [ax] = graph1D_jd(rho,gradpthermalsgradrho(:,1),0,0,'\rho','\nablap\cdot\psi/|\nabla\rho|',strrep(equil.id,'_','-'),NaN,NaN,NaN,'-','none','r',2);drawnow
0403 end
0404 [ax] = graph1D_jd(equil.eqdsk.prho,equil.eqdsk.gradpkinsgradrho,0,0,'\rho','\nablap\cdot\psi/|\nabla\rho|',strrep(equil.id,'_','-'),leg,NaN,NaN,'-','none','g',1);
0405 else
0406 [ax] = graph1D_jd(rho,gradpmagsgradrho,0,0,'\rho','\nablap\cdot\psi/|\nabla\rho|',strrep(equil.id,'_','-'),NaN,NaN,NaN,'-','none','b',2);drawnow
0407 if ~isnan(jparz_ohm_fsav),
0408 [ax] = graph1D_jd(rho,gradpthermalsgradrho(:,1),0,0,'\rho','\nablap\cdot\psi/|\nabla\rho|',strrep(equil.id,'_','-'),leg,NaN,NaN,'-','none','r',2);drawnow
0409 end
0410 end
0411 print_jd(2,['Pressure_gradient_',equil.id]);
0412
0413 figure('Name','Axial current density'),clf
0414 [ax] = graph1D_jd(rho,jz/1e6,0,0,'\rho','j_{||z} (MA/m^{2})',[strrep(equil.id,'_','-'),', \theta\in[0:\pi/',int2str((mtheta-1)/2),':2\pi], I_{p} = ',num2str(Ip_z/1e6),' (MA), li = ',num2str(li)],NaN,NaN,NaN,'-','none','r',2);drawnow
0415 [ax] = graph1D_jd(rho,jz_fsav/1e6,0,0,'\rho','j_{||z} (MA/m^{2})',[strrep(equil.id,'_','-'),', \theta\in[0:\pi/',int2str((mtheta-1)/2),':2\pi]'],NaN,NaN,NaN,'-','none','b',1);drawnow
0416 if isfield(equil,'eqdsk'),
0417
0418 end
0419
0420 leg = {'j_{z}','<j_{z}>','<j_{z}>_{(eqdsk)}','Location','NorthEast'};
0421 axf = axis(ax);
0422 text(0.1,axf(3)+(axf(4)-axf(3))/10, ['I_p = ',num2str(Ip_z/1e6),' (MA)'], 'Color', 'b','FontSize',18);
0423 text(0.1,axf(3)+2*(axf(4)-axf(3))/10, ['l_i = ',num2str(li),], 'Color', 'b','FontSize',18);
0424 print_jd(2,['Axial_current_density_',equil.id]);
0425
0426 figure('Name','Axial component of the parallel current density'),clf
0427 [ax] = graph1D_jd(rho,jparz/1e6,0,0,'\rho','j_{||z} (MA/m^{2})',[strrep(equil.id,'_','-'),', \theta\in[0:\pi/',int2str((mtheta-1)/2),':2\pi], I_{p} = ',num2str(Ip_z/1e6),' (MA), li = ',num2str(li)],NaN,NaN,NaN,'-','none','r',2);drawnow
0428 [ax] = graph1D_jd(rho,jparz_fsav/1e6,0,0,'\rho','j_{||z} (MA/m^{2})',[strrep(equil.id,'_','-'),', \theta\in[0:\pi/',int2str((mtheta-1)/2),':2\pi]'],NaN,NaN,NaN,'-','none','b',1);drawnow
0429 if isfield(equil,'eqdsk'),
0430
0431 end
0432
0433 leg = {'j_{||z}','<j_{||z}>','Location','NorthEast'};
0434 axf = axis(ax);
0435 text(0.1,axf(3)+(axf(4)-axf(3))/10, ['I_p = ',num2str(Ip_parperpz/1e6),' (MA)'], 'Color', 'b','FontSize',18);
0436 text(0.1,axf(3)+2*(axf(4)-axf(3))/10, ['l_i = ',num2str(li),], 'Color', 'b','FontSize',18);
0437 print_jd(2,['Parallel_z_current_density_',equil.id]);
0438
0439 figure('Name','Axial component of the perpendicular current density'),clf
0440 [ax] = graph1D_jd(rho,jperpz/1e6,0,0,'\rho','j_{||z} (MA/m^{2})',[strrep(equil.id,'_','-'),', \theta\in[0:\pi/',int2str((mtheta-1)/2),':2\pi], I_{p} = ',num2str(Ip_z/1e6),' (MA), li = ',num2str(li)],NaN,NaN,NaN,'-','none','r',2);drawnow
0441 [ax] = graph1D_jd(rho,jperpz_fsav/1e6,0,0,'\rho','j_{||z} (MA/m^{2})',[strrep(equil.id,'_','-'),', \theta\in[0:\pi/',int2str((mtheta-1)/2),':2\pi]'],NaN,NaN,NaN,'-','none','b',1);drawnow
0442
0443 leg = {'j_{perpz}','<j_{perpz}>','Location','NorthEast'};
0444 axf = axis(ax);
0445 text(0.1,axf(3)+(axf(4)-axf(3))/10, ['I_p = ',num2str(Ip_parperpz/1e6),' (MA)'], 'Color', 'b','FontSize',18);
0446 text(0.1,axf(3)+2*(axf(4)-axf(3))/10, ['l_i = ',num2str(li),], 'Color', 'b','FontSize',18);
0447 print_jd(2,['Perpendicular_z_current_density_',equil.id]);
0448
0449 if isfield(equil,'eqdsk'),
0450
0451 figure('Name','Safety factor'),clf
0452 leg = {'Magnetic','EQDSK','Location','NorthEast'};
0453 if isinf(equilDKE.Rp),
0454 [ax] = graph1D_jd(equil.eqdsk.prho,equil.eqdsk.pq*equilDKE.Rp/equilDKE.ap,0,0,'\rho','q(\rho)*R_p/a_p',strrep(equil.id,'_','-'),NaN,NaN,NaN,'-','none','r',2);drawnow
0455 [ax] = graph1D_jd(equil.eqdsk.prho(2:end),xq,0,0,'\rho','q(\rho)*R_p/a_p',strrep(equil.id,'_','-'),leg,NaN,NaN,'-','none','b',1);drawnow
0456 else
0457 [ax] = graph1D_jd(equil.eqdsk.prho,equil.eqdsk.pq,0,0,'\rho','q(\rho)',strrep(equil.id,'_','-'),NaN,NaN,NaN,'-','none','r',2);drawnow
0458 [ax] = graph1D_jd(equil.eqdsk.prho(2:end),xq*equilDKE.ap/equilDKE.Rp,0,0,'\rho','q(\rho)',strrep(equil.id,'_','-'),leg,NaN,NaN,'-','none','b',1);drawnow
0459 end
0460 axf = axis(ax);
0461 text(0.1,axf(3)+9*(axf(4)-axf(3))/10, ['I_p = ',num2str(Ip_z/1e6),' (MA)'], 'Color', 'b','FontSize',18);
0462 text(0.1,axf(3)+8*(axf(4)-axf(3))/10, ['l_i = ',num2str(li),], 'Color', 'b','FontSize',18);
0463 print_jd(2,['Safety_factor_',equil.id]);
0464
0465 figure('Name','Current flux function'),clf
0466
0467 if ~isinf(equilDKE.Rp),
0468 [ax] = graph1D_jd(rho,pF_eqdsk(:,1),0,0,'\rho','F(\psi)=RB_{\phi} [mT]',strrep(equil.id,'_','-'),NaN,NaN,NaN,'-','none','r',2);drawnow
0469 [ax] = graph1D_jd(rho,F(:,1),0,0,'\rho','F(\psi)=RB_{\phi} [mT]',strrep(equil.id,'_','-'),leg,NaN,NaN,'-','none','b',1);drawnow
0470 leg = {'EQDSK','Magnetic','Location','NorthEast'};
0471 [ax] = graph1D_jd(rho,pF_eqdsk(:,2:end),0,0,'\rho','F(\psi)=RB_{\phi} [mT]',strrep(equil.id,'_','-'),NaN,NaN,NaN,'-','none','r',2);drawnow
0472 [ax] = graph1D_jd(rho,F(:,2:end),0,0,'\rho','F(\psi)=RB_{\phi} [mT]',strrep(equil.id,'_','-'),leg,NaN,NaN,'-','none','b',1);drawnow
0473 end
0474 axf = axis(ax);
0475 text(0.1,axf(3)+9*(axf(4)-axf(3))/10, ['I_p = ',num2str(Ip_z/1e6),' (MA)'], 'Color', 'b','FontSize',18);
0476 text(0.1,axf(3)+8*(axf(4)-axf(3))/10, ['l_i = ',num2str(li),], 'Color', 'b','FontSize',18);
0477 print_jd(2,['Current_flux_function_',equil.id]);
0478
0479 figure('Name','Flux derivative of the current flux function'),clf
0480 leg = {'Magnetic','EQDSK','Location','NorthEast'};
0481
0482 if ~isinf(equilDKE.Rp),
0483 [ax] = graph1D_jd(rho,pFdFdpsi_eqdsk,0,0,'\rho','q(\rho)',strrep(equil.id,'_','-'),NaN,NaN,NaN,'-','none','r',2);drawnow
0484 [ax] = graph1D_jd(rho,dFdpsi,0,0,'\rho','F(\psi)',strrep(equil.id,'_','-'),leg,NaN,NaN,'-','none','b',1);drawnow
0485 end
0486 axf = axis(ax);
0487 text(0.1,axf(3)+9*(axf(4)-axf(3))/10, ['I_p = ',num2str(Ip_z/1e6),' (MA)'], 'Color', 'b','FontSize',18);
0488 text(0.1,axf(3)+8*(axf(4)-axf(3))/10, ['l_i = ',num2str(li),], 'Color', 'b','FontSize',18);
0489 print_jd(2,['Flux_derivative _urrent_flux_function_',equil.id]);
0490
0491
0492 end
0493
0494 if ~isnan(jparz_ohm_fsav),
0495
0496 figure('Name','Current density and li'),clf
0497 [ax] = graph1D_jd(rho,jparz_fsav/1e6/max(abs(jparz_fsav/1e6)),0,0,'\rho','j_{||z}/max(j_{||z})',[strrep(equil.id,'_','-'),', \theta\in[0:\pi/',int2str((mtheta-1)/2),':2\pi]'],NaN,NaN,NaN,'-','none','b',2);drawnow
0498 [ax] = graph1D_jd(rho,jparz_ohm_fsav/1e6,0,0,'\rho','j_{||z}/max(j_{||z})',[strrep(equil.id,'_','-'),', \theta\in[0:\pi/',int2str((mtheta-1)/2),':2\pi]'],leg,NaN,NaN,'-','none','g',1);drawnow
0499
0500 leg = {['<j_{||z}>, l_i=',num2str(li)],'<j_{||z}>_{Ohm}','Location','NorthEast'};
0501
0502 for j = 1:6,
0503 [ax] = graph1D_jd(rho,jparz_fsav_cyl(:,j)/max(abs(jparz_fsav_cyl(:,j))),0,0,'\rho','j_{||z}/max(j_{||z})',[strrep(equil.id,'_','-'),', \theta\in[0:\pi/',int2str((mtheta-1)/2),':2\pi], I_{p} = ',num2str(Ip/1e6),' (MA), li = ',num2str(li)],leg,NaN,NaN,'--','none','r',1);drawnow
0504 end
0505
0506 text(0.02,0.52, '<j_{||z}>_{circ.} ~ (1-\rho^2)^n', 'Color', 'r','FontSize',18);
0507 text(0.02,0.44, ['n = 5, l_i = ',num2str(1.80)], 'Color', 'r','FontSize',18);
0508 text(0.02,0.36, ['n = 4, l_i = ',num2str(1.67)], 'Color', 'r','FontSize',18);
0509 text(0.02,0.28, ['n = 3, l_i = ',num2str(1.45)], 'Color', 'r','FontSize',18);
0510 text(0.02,0.20, ['n = 2, l_i = ',num2str(1.21)], 'Color', 'r','FontSize',18);
0511 text(0.02,0.12, ['n = 1, l_i = ',num2str(0.91)], 'Color', 'r','FontSize',18);
0512 text(0.02,0.04, ['n = 0, l_i = ',num2str(0.5)], 'Color', 'r','FontSize',18);
0513
0514 print_jd(2,['Current_density_li_',equil.id]);
0515 end
0516
0517 disp(['Normalized surface Sn = ',num2str(Surf/Surf_cyl)])
0518
0519 end
0520
0521 equil.ip = Ip_z/1e6;
0522 equil.li = li;
0523 equil.shafranov_shift = pshift(end);
0524 equil.elongation = pelong(end);
0525 equil.triangularity = ptrian(end);
0526 equil.Te_vol = Te_vol;
0527 equil.ne_vol = ne_vol;
0528 equil.zTi_vol = zTi_vol;
0529 equil.zni_vol = zni_vol;
0530
0531 equil.consistency.rho = rho;
0532 equil.consistency.jparz_fsav = jparz_fsav;
0533 equil.consistency.jparz_fsav_cyl = jparz_fsav_cyl;
0534 equil.consistency.Te = Te;
0535 equil.consistency.ne = ne;
0536
0537 if isfield(equil,'eqdsk')
0538 equil.consistency.rhoP = rhoP;
0539 equil.consistency.rhoT = rhoT;
0540 equil.consistency.rhoV = rhoV;
0541
0542 equil.consistency.q = q;
0543 end
0544
0545