Description of eqdsk

0001 function [x,y,xypsin,psia,rpsin,rF,rp,rq,Rp,Zp,Ip,Rv,Zv,Bv,rFdFdpsi,rdpdpsi] = eqdsk_jd(eqdsk,display_mode,p_opt,dispname,savename,Ns)
0002 % EFIT - function that reads the equilibrium parameters from eqdsk files by EFIT
0003 %
0004 % This function reads the equilibrium parameters from eqdsk files created
0005 % by EFIT and saved in the formatted text file "filename"
0006 %
0007 % INPUT
0008 %
0009 %    - eqdsk: file that contains the EQDSK structure
0010 %    - display_mode: display option (0,2)
0011 %       (default display_mode = 0, no display)
0012 %    - p_opt: printing option
0013 %       (default p_opt = 2, print mutiple formats)
0014 %    - dispname: name displayed in the figures
0015 %       (default dispname = '')
0016 %    - savename: name used for saving the the figures
0017 %       (default savename = 'equil')
0018 %    - Ns: magnetic axis recalculation supergrid factor
0019 %       (default Ns = 101;
0020 %
0021 % OUTPUT:
0022 %
0023 %    - x: R - Rp values on the nnr * nnz equispaced grid
0024 %    - y: Z - Zp values on the nnr * nnz equispaced grid
0025 %    - xypsin: normalized psi values on the nnr * nnz equispaced grid
0026 %    - psia: poloidal flux function value at the last closed flux surface (touching the limiter or the separatrix)
0027 %    - rpsin: normalized poloidal magnetic flux grid
0028 %    - rF: R*BPHI at nnr equispaced points in psi from psimag to psilim.
0029 %    - rp: pressure values at nnr equispaced points in psi from psimag to psilim.
0030 %    - rq: safety factor q on the equispaced psi grid
0031 %    - Rp: major radius of magnetic axis
0032 %    - Zp: vertical height of magnetic axis
0033 %    - Ip: toroidal lasma current
0034 %    - Rv: nominal major radius of the torus
0035 %    - Zv: vertical shift of the rectangular box midplane
0036 %    - Bv: vacuum toroidal magnetic field at radmaj
0037 %    - rFdFdpsi: ffpsiar_prime, [i.e., d f/d (psi)] at nnr equispaced points in psi from psimag to psilim.
0038 %    - rdpdpsi: p_prime [i.e., d prar/d (psi)] at nnr equispaced points in psi from psimag to psilim.
0039 %
0040 % WARNING: EFIT sign treatment
0041 %
0042 % By convention, the EFIT calculations always give psia > 0 regardless of the direction of the
0043 % toroidal current. The signs of Bv and Ip are given with respect to phi_eqdsk which is such that
0044 % (R,phi_eqdsk,Z) is right handed. The sign of phi_dke is opposite such that (R,Z,phi_dke) and
0045 % (r,theta,phi_dke) be both right handed. The outputs of eqdsk_jd must respect the DKE convention.
0046 % In that case the signs of Ip and of psia must be identical.
0047 %
0048 % by J. Decker (joan.decker@alum.mit.edu) 10/28/04
0049 %
0050 %
0051 if nargin < 6
0052     Ns = 101;%magnetic axis recalculation supergrid factor
0053 end
0054 %
0055 if nargin < 5
0056     savename = 'equil';
0057 end
0058 %
0059 if nargin < 4
0060     dispname = '';
0061 end
0062 %
0063 if nargin < 3
0064     p_opt = 2;
0065 end
0066 %
0067 if nargin < 2
0068     display_mode = 0;
0069 end
0070 %
0071 if ~isstruct(eqdsk),
0072     %
0073     if exist(eqdsk,'file'),
0074         eqdsk = fileread(eqdsk);
0075     end
0076     %
0077     istr = find(double(eqdsk) == 10,1,'first');
0078     str_intro = eqdsk(1:istr-1);
0079     eqdsk = eqdsk(istr+1 : end);
0080     sc = 0;
0081     nn = [];
0082     while length(str_intro) > 0 & length(nn) < 3,
0083         if str_intro(end) ~= ' ' 
0084             if sc == 0,
0085                 sc = 1;
0086                 str = str_intro(end);
0087             else
0088                 str = [str_intro(end),str];
0089             end
0090         else
0091             if sc == 1,
0092                 nn = [nn,str2num(str)];
0093                 sc = 0;
0094             end           
0095         end
0096         str_intro(end) = '';
0097     end
0098     %
0099     if nn(3) < 10,
0100         nnx = nn(2);
0101         nny = nn(1);
0102         nnr = nn(2);
0103     else
0104         nnx = nn(3);
0105         nny = nn(2);
0106         if nn(1) == 0,
0107             nnr = nn(2);%CQL3D (COMPASS)
0108         else
0109             nnr = nn(1);
0110         end
0111     end
0112     %
0113     [data,~,~,ind] = sscanf(eqdsk,'%e',5);
0114     xbox = data(1);% horizontal full-width of the rectangle on which poloidal flux values are given
0115     ybox = data(2);% vertical full-width of the rectangle on which poloidal flux values are given
0116     Rv = data(3);% nominal major radius of the torus
0117     Rmin = data(4);% major radius of the inner edge of the box of rectangular grid
0118     Zv = data(5);% vertical shift of the rectangular box midplane
0119     %
0120     eqdsk = eqdsk(ind:end);
0121     [data,~,~,ind] = sscanf(eqdsk,'%e',5);
0122     Rp = data(1);% major radius of magnetic axis
0123     Zp = data(2);% vertical height of magnetic axis
0124     psip = data(3);% poloidal flux function value at the magnetic axis
0125     psia = data(4);% poloidal flux function value at the last closed flux surface (touching the limiter or the separatrix)
0126     Bv = data(5);% vacuum toroidal magnetic field at radmaj
0127     %
0128     eqdsk = eqdsk(ind:end);
0129     [data,~,~,ind] = sscanf(eqdsk,'%e',5);
0130     Ip = data(1);% toroidal current
0131     psimx1 = data(2);% not used
0132     psimx2 = data(3);% not used
0133     xax1 = data(4);% not used
0134     xax2 = data(5);% not used
0135     %
0136     eqdsk = eqdsk(ind:end);
0137     [data,~,~,ind] = sscanf(eqdsk,'%e',5);
0138     zax1 = data(1);% not used
0139     zax2 = data(2);% not used
0140     psisep = data(3);% not used
0141     xsep = data(4);% not used
0142     zsep = data(5);% not used
0143     %
0144     eqdsk = eqdsk(ind:end);
0145     [rF,~,~,ind] = sscanf(eqdsk,'%e',nnr);% R*BPHI  at nnr equispaced points in psi from psimag to psilim.
0146     eqdsk = eqdsk(ind:end);
0147     [rp,~,~,ind] = sscanf(eqdsk,'%e',nnr);% pressure values at the same points
0148     eqdsk = eqdsk(ind:end);
0149     [rFdFdpsi,~,~,ind] = sscanf(eqdsk,'%e',nnr);% ffpsiar_prime, [i.e., d f/d (psi)] at the same points
0150     eqdsk = eqdsk(ind:end);
0151     [rdpdpsi,~,~,ind] = sscanf(eqdsk,'%e',nnr);% p_prime [i.e., d prar/d (psi)] at the same points
0152     %
0153     rF = rF.';
0154     rp = rp.';
0155     rFdFdpsi = rFdFdpsi.';
0156     rdpdpsi = rdpdpsi.';
0157     %
0158     eqdsk = eqdsk(ind:end);
0159     [xypsi,~,~,ind] = sscanf(eqdsk,'%e',[nnx,nny]);% psi values on the nnr * nnz equispaced grid
0160     %
0161     eqdsk = eqdsk(ind:end);
0162     [rq,~,~,ind] = sscanf(eqdsk,'%e',nnr); % safety factor q on the equispaced psi grid
0163     rq = rq.';
0164     %
0165     eqdsk = eqdsk(ind:end);
0166     [ncoil,~,~,ind] = sscanf(eqdsk,'%d',5);
0167     ncoil = ncoil(ncoil>0).';% not used
0168     eqdsk = eqdsk(ind:end);
0169     [ccoil,~,~,ind] = sscanf(eqdsk,'%e',[ncoil(1),length(ncoil)]);% not used
0170 else
0171     nnx = eqdsk.nnx;
0172     nny = eqdsk.nny;
0173     nnr = eqdsk.nnr;
0174     xbox = eqdsk.xbox;% horizontal full-width of the rectangle on which poloidal flux values are given
0175     ybox = eqdsk.ybox;% vertical full-width of the rectangle on which poloidal flux values are given
0176     Rv = eqdsk.Rv;% nominal major radius of the torus
0177     Rmin = eqdsk.Rmin;% major radius of the inner edge of the box of rectangular grid
0178     Zv = eqdsk.Zv;% vertical shift of the rectangular box midplane
0179     %
0180     Rp = eqdsk.Rp;% major radius of magnetic axis
0181     Zp = eqdsk.Zp;% vertical height of magnetic axis
0182     psip = eqdsk.psip;% poloidal flux function value at the magnetic axis
0183     psia = eqdsk.psia;% poloidal flux function value at the last closed flux surface (touching the limiter or the separatrix)
0184     Bv = eqdsk.Bv;% vacuum toroidal magnetic field at radmaj
0185     %
0186     Ip = eqdsk.Ip;% toroidal current
0187     %
0188     rF = eqdsk.rF;% R*BPHI  at nnr equispaced points in psi from psimag to psilim.
0189     rp = eqdsk.rp;% pressure values at the same points
0190     rFdFdpsi = eqdsk.rFdFdpsi;% ffpsiar_prime, [i.e., d f/d (psi)] at the same points
0191     rdpdpsi = eqdsk.rdpdpsi;% p_prime [i.e., d prar/d (psi)] at the same points
0192     %
0193     rF = rF.';
0194     rp = rp.';
0195     rFdFdpsi = rFdFdpsi.';
0196     rdpdpsi = rdpdpsi.';
0197     %
0198     xypsi = eqdsk.xypsi;% psi values on the nnr * nnz equispaced grid
0199     %
0200     rq = eqdsk.rq; % safety factor q on the equispaced psi grid
0201     rq = rq.';
0202 end
0203 %
0204 R = Rmin + linspace(0,xbox,nnx);
0205 Z = Zv + linspace(-ybox/2,ybox/2,nny);
0206 %
0207 % To compensate irregularities in the EFIt output
0208 % Rp, Zp and psip are recalculated in accordance
0209 % with xypsi data
0210 %
0211 psip_old = psip;
0212 Rp_old = Rp;
0213 Zp_old = Zp;
0214 %
0215 if psip < psia, %IPHI is positive, PHI clockwise from top
0216     sigma = 1;
0217 else %IPHI is negative, PHI clockwise from top
0218     sigma = -1;
0219 end
0220 %
0221 % create mask to search for minimum... within the plasma
0222 %
0223 xR = repmat(R.',[1,nny]);
0224 xZ = repmat(Z,[nnx,1]);
0225 xd = sqrt((xR - Rp).^2 + (xZ - Zp).^2);
0226 xmask_out = (sigma*xypsi > sigma*psia);% these points are outside the plasma
0227 dmin = min(xd(xmask_out));
0228 xmask_excl = (xd > dmin);
0229 xypsi_excl = xypsi;
0230 xypsi_excl(xmask_excl) = sigma*Inf;
0231 %
0232 [iR,iZ] = find(sigma*xypsi_excl == min(min(sigma*xypsi_excl)));
0233 %
0234 if length(iR) ~= 1,
0235     error('unable to locate minimum in xypsi data');
0236 end
0237 %
0238 ss = linspace(0,1,Ns);
0239 sR = R(iR-1) + (R(iR+1)-R(iR-1))*ss;
0240 sZ = Z(iZ-1) + (Z(iZ+1)-Z(iZ-1))*ss;
0241 %
0242 sspsi = interp2(R,Z,xypsi.',sR.',sZ,'spline').';
0243 psip = sigma*min(min(sigma*sspsi));
0244 [isR,isZ] = find(sspsi == psip);
0245 %
0246 if length(isR) ~= 1,
0247     error('unable to locate minimum in xypsi data');
0248 end
0249 %
0250 Rp = sR(isR);
0251 Zp = sZ(isZ);
0252 %
0253 if display_mode >= 1,
0254     disp(['Old value of psip: ',num2str(psip_old),', new value: ',num2str(psip)])
0255     disp(['Old value of Rp: ',num2str(Rp_old),', new value: ',num2str(Rp)])
0256     disp(['Old value of Zp: ',num2str(Zp_old),', new value: ',num2str(Zp)])
0257 end
0258 %
0259 % iRp = find(abs(R - Rp) == min(abs(R - Rp)));
0260 % if abs(R(iRp) - Rp)/Rp < 1/(100*nnr) %one R point is very close to Rp
0261 %     R = R + Rp - R(iRp);
0262 %     psimag = min(min(epsi));
0263 % else
0264 %     warning('No R = Rp point');
0265 % end
0266 %
0267 %
0268 x = R - Rp;
0269 y = Z - Zp;
0270 %
0271 xypsi = xypsi - psip;
0272 psia = (psia - psip);
0273 %
0274 xypsin = xypsi/psia;
0275 rpsin = linspace(0,1,nnr);
0276 %
0277 % EFIT sign treatment
0278 %
0279 % By convention, the EFIT calculations always give psia > 0 regardless of the direction of the
0280 % toroidal current. The signs of Bv and Ip are given with respect to phi_eqdsk which is such that
0281 % (R,phi_eqdsk,Z) is right handed. The sign of phi_dke is opposite such that (R,Z,phi_dke) and
0282 % (r,theta,phi_dke) be both right handed. The outputs of eqdsk_jd must respect the DKE convention.
0283 % In that case the signs of Ip and of psia must be identical.
0284 %
0285 Ip = - Ip;
0286 Bv = - Bv;
0287 rF = - rF;
0288 rFdFdpsi = - rFdFdpsi;
0289 %
0290 if sign(Ip) ~= sign(psia),
0291     psia = - psia;
0292 end
0293 %
0294 if sign(Bv) ~= sign(rF(1)),
0295     rF = - rF;
0296     rFdFdpsi = - rFdFdpsi;
0297 end
0298 %
0299 if display_mode >= 2,
0300     %
0301     style = '-';
0302     marker = '+';
0303     color = 'r';
0304     width = 2;
0305     siz = 20;
0306     red = 0.9;
0307     lspace = 0.7;
0308     bspace = 0.6;
0309     tit = strrep(dispname,'_','-');
0310     %
0311     figure(1),clf
0312     %
0313     xlab = 'x (m)';
0314     ylab = 'y (m)';
0315     xlim = [min(x),max(x)];
0316     ylim = [min(y),max(y)];
0317     %
0318     graph2D_jd(x,y,xypsin,xlab,ylab,tit,xlim,ylim,0,NaN,[0:0.05:1],0,style,color,width,siz,red);
0319     axis('equal');
0320     axis([xlim,ylim]);
0321     %set(gca,'XTick',[0:1:2],'YTick',[-2:2]);
0322     print_jd(p_opt,[savename,'_psi_xy']);
0323     %
0324     figure(2),clf
0325     %
0326     xlab = '\psi_n^{1/2}';
0327     ylab = 'F(\psi)=RB_{\phi} (mT)';
0328     xlim = [0,1];
0329     ylim = NaN;
0330     %
0331     graph1D_jd(sqrt(rpsin),rF,0,0,xlab,ylab,tit,NaN,xlim,ylim,style,marker,color,width,siz,gca,red,lspace,bspace);
0332     set(gca,'XTick',[0:0.2:1]);
0333     print_jd(p_opt,[savename,'_F_psi']);
0334     %
0335     figure(3),clf
0336     %
0337     ylab = 'p(\psi) (hPa)';
0338     ylim = NaN;
0339     %
0340     graph1D_jd(sqrt(rpsin),rp/1e2,0,0,xlab,ylab,tit,NaN,xlim,ylim,style,marker,color,width,siz,gca,red,lspace,bspace);
0341     set(gca,'XTick',[0:0.2:1]);
0342     print_jd(p_opt,[savename,'_p_psi']);
0343     %
0344     figure(4),clf
0345     %
0346     ylab = 'FdF/d\psi (m^{-1})';
0347     ylim = NaN;
0348     %
0349     graph1D_jd(sqrt(rpsin),rFdFdpsi,0,0,xlab,ylab,tit,NaN,xlim,ylim,style,marker,color,width,siz,gca,red,lspace,bspace);
0350     set(gca,'XTick',[0:0.2:1]);
0351     print_jd(p_opt,[savename,'_FdF_psi']);
0352     %
0353     figure(5),clf
0354     %
0355     ylab = 'dp/d\psi (hPaB^{-1}m^{-2})';
0356     ylim = NaN;
0357     %
0358     graph1D_jd(sqrt(rpsin),rdpdpsi/1e5,0,0,xlab,ylab,tit,NaN,xlim,ylim,style,marker,color,width,siz,gca,red,lspace,bspace);
0359     set(gca,'XTick',[0:0.2:1]);
0360     print_jd(p_opt,[savename,'_dp_psi']);
0361     %
0362     figure(6),clf
0363     %
0364     ylab = 'q(\psi)';
0365     ylim = NaN;
0366     %
0367     graph1D_jd(sqrt(rpsin),rq,0,0,xlab,ylab,tit,NaN,xlim,ylim,style,marker,color,width,siz,gca,red,lspace,bspace);
0368     set(gca,'XTick',[0:0.2:1]);
0369     print_jd(p_opt,[savename,'_q_psi']);
0370     %
0371     figure(7),clf
0372     %
0373     rho_list = [0:0.05:0.999];
0374     iy0 = find(abs(y) == min(abs(y)));
0375     ix0 = min(find(x >= 0));
0376     ap = interp1(xypsin(ix0:nnx,iy0(1)),x(ix0:nnx),1,'spline');
0377     rpsi = interp1(x(ix0:nnx)/ap,xypsin(ix0:nnx,iy0(1)),rho_list,'spline');
0378     %
0379     xlab = 'x (m)';
0380     ylab = 'y (m)';
0381     xlim = [min(x),max(x)];
0382     ylim = [min(y),max(y)];
0383     tit = strrep(dispname,'_','-');
0384     %
0385     graph2D_jd(x,y,xypsin,xlab,ylab,tit,xlim,ylim,0,NaN,rpsi,0,style,color,width,siz,red);
0386     axis('equal');
0387     axis([xlim,ylim]);
0388     %set(gca,'XTick',[0:1:2],'YTick',[-2:2]);
0389     print_jd(p_opt,[savename,'_rho_xy']);
0390  end
0391 
0392 
0393 
0394 
0395 
0396 
0397 
0398
eqdsk_jd

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SOURCE CODE
