equil_EFIT_jd

 This function reads the equilibrium parameters from eqdsk files created
 by EFIT and saved in the formatted text file "filename.efit", as well as
 the temperature and density parameters from the (fully numerical) text
 file "filename". It creates a DKE equilibrium file "EQUIL_filename.mat"

    INPUTS:

        - filename: equilibrium data and save name
        - type: profile data type:(0) for ascii; (1) for .mat
        - display_mode: option for figure display
        - p_opt: option for printing

 by J. Decker (jodecke@alum.mit.edu) 11/01/04

0001 function equil = equil_EFIT_jd(filename,type,display_mode,p_opt,dispname,savename,Zeff)
0002 %
0003 % This function reads the equilibrium parameters from eqdsk files created
0004 % by EFIT and saved in the formatted text file "filename.efit", as well as
0005 % the temperature and density parameters from the (fully numerical) text
0006 % file "filename". It creates a DKE equilibrium file "EQUIL_filename.mat"
0007 %
0008 %    INPUTS:
0009 %
0010 %        - filename: equilibrium data and save name
0011 %        - type: profile data type:(0) for ascii; (1) for .mat
0012 %        - display_mode: option for figure display
0013 %        - p_opt: option for printing
0014 %
0015 % by J. Decker (jodecke@alum.mit.edu) 11/01/04
0016 %
0017 if nargin < 7,% & type == 3
0018     Zeff = NaN;
0019 end
0020 if nargin < 6,
0021     savename = filename;
0022 end
0023 if nargin < 5,
0024     dispname = filename;
0025 end
0026 method = 'spline';%'linear';
0027 nt = 65;
0028 %
0029 % Reads EFIT eqdsk data
0030 %
0031 efitname = [filename,'.efit'];
0032 [x,y,xypsin,psia,psin,pF,pp,pq,Rp,Zp,Ip,Rv,Zv,Bv,pFdFdpsi,pdpdpsi] = eqdsk_jd(efitname,display_mode,p_opt,dispname,savename);
0033 %
0034 % Reads EFIT numerical text T and n data
0035 %
0036 profname = [filename,'.prof'];
0037 if exist(profname,'file')
0038     [yrho,yspsin,yspsinT,yTe,yne,yzTi,yzni,zZi,zmi] = Teneprof_jd(profname,type,display_mode,p_opt,Zeff,dispname,savename);
0039     flagprof = 1;
0040 else
0041     flagprof = 0;
0042 end
0043 %
0044 [ptx,pty,ptdpsindx,ptdpsindy,psin,theta] = rz2pt_jd(x,y,xypsin,psin,nt,method,display_mode,p_opt,savename);
0045 %
0046 np = length(psin);
0047 ap = ptx(np,1);%Plasma minor radius on the LFS midplane
0048 psi_apRp = psin*psia*ap/Rp;%The poloidal flux coordinate is normalized by the aspect ratio
0049 ptR = Rp + ptx;
0050 %
0051 % Note that in DKE we have B = F.gradphi + gradphi x gradpsi.
0052 %
0053 ptBx = -ptdpsindy*psia./ptR;
0054 ptBy = ptdpsindx*psia./ptR;
0055 %
0056 ptF = repmat(pF.',[1,nt]);
0057 ptBPHI = ptF./ptR;
0058 %
0059 prho = ptx(:,1).'/ap;
0060 %
0061 pq_h = (pq(1:end-1) + pq(2:end))/2;
0062 pdpsin = diff(psin);
0063 ppsiT = [0,cumsum(pq_h.*pdpsin)];
0064 pspsinT = sqrt(ppsiT/ppsiT(end));
0065 %
0066 pspsin = sqrt(psin);
0067 %
0068 if flagprof && ~isempty(yrho),% profiles given on rho grid
0069     pTe = interp1(yrho,yTe,prho,method);
0070     pne = interp1(yrho,yne,prho,method);
0071     pzTi = interp1(yrho,yzTi.',prho,method).';
0072     pzni = interp1(yrho,yzni.',prho,method).';
0073     %
0074 elseif flagprof && ~isempty(yspsin),% profiles given on poloidal flux based grid
0075     pTe = interp1(yspsin,yTe,pspsin,method);
0076     pne = interp1(yspsin,yne,pspsin,method);
0077     pzTi = interp1(yspsin,yzTi.',pspsin,method).';
0078     pzni = interp1(yspsin,yzni.',pspsin,method).'; 
0079     %
0080 elseif flagprof &&~isempty(yspsinT),% profiles given on toroidal flux based grid
0081     pTe = interp1(yspsinT,yTe,pspsinT,method);
0082     pne = interp1(yspsinT,yne,pspsinT,method);
0083     pzTi = interp1(yspsinT,yzTi.',pspsinT,method).';
0084     pzni = interp1(yspsinT,yzni.',pspsinT,method).';    
0085     %
0086 end
0087 %
0088 if display_mode == 2, 
0089 %
0090     rx = [-fliplr(prho(2:np)),prho];
0091     rBx = [fliplr(ptBx(2:np,(nt+1)/2)'),ptBx(:,1)'];
0092     rBy = [fliplr(ptBy(2:np,(nt+1)/2)'),ptBy(:,1)'];
0093     rBPHI = [fliplr(ptBPHI(2:np,(nt+1)/2)'),ptBPHI(:,1)'];
0094     rBP = sqrt(rBx.^2 + rBy.^2);
0095     rB = sqrt(rBP.^2 + rBPHI.^2);
0096 %
0097     style = '-';
0098     marker = 'none';
0099     color1 = 'r';
0100     color2 = 'b';
0101     color3 = 'k';
0102     width = 2;
0103     siz = 20;
0104     red = 0.9;
0105     lspace = 0.7;
0106     bspace = 0.6;
0107     tit = [dispname];
0108     %
0109     figure(12),clf
0110     %
0111     xlab = 'r/a';
0112     ylab = 'B(T)';
0113     xlim = [-1,1];
0114     ylim = NaN;
0115     leg = ['B_T';'B_P';'B  '];
0116 %
0117     graph1D_jd(rx,abs(rBPHI),0,0,'','','',NaN,xlim,ylim,style,marker,color1,width,siz,gca,red,lspace,bspace);
0118     graph1D_jd(rx,rBP,0,0,'','','',NaN,xlim,ylim,style,marker,color2,width);
0119     graph1D_jd(rx,rB,0,0,xlab,ylab,tit,leg,xlim,ylim,style,marker,color3,width);
0120     print_jd(p_opt,[savename,'_Bfield']);
0121 end
0122 %
0123 ptx(:,65) = ptx(:,1);
0124 pty(:,65) = pty(:,1);
0125 ptBx(:,65) = ptBx(:,1);
0126 ptBy(:,65) = ptBy(:,1);
0127 ptBPHI(:,65) = ptBPHI(:,1);
0128 %
0129 % equil structure
0130 %
0131 equil.id = savename;
0132 %
0133 if flagprof,
0134     equil.zZi = zZi;
0135     equil.zmi = zmi;
0136     equil.pTe = pTe;
0137     equil.pne = pne;
0138     equil.pzTi = pzTi;
0139     equil.pzni = pzni;
0140 end
0141 %
0142 equil.psi_apRp = psi_apRp;
0143 equil.theta = theta;
0144 equil.ptx = ptx;
0145 equil.pty = pty;
0146 equil.ptBx = ptBx;
0147 equil.ptBy = ptBy;
0148 equil.ptBPHI = ptBPHI;
0149 equil.Rp = Rp;
0150 equil.Zp = Zp;
0151 %
0152 equil_magnetic.eqdsk.x = x;
0153 equil_magnetic.eqdsk.y = y;
0154 equil_magnetic.eqdsk.xypsin = xypsin;
0155 equil_magnetic.eqdsk.psia = psia;
0156 equil_magnetic.eqdsk.psin = psin;
0157 equil_magnetic.eqdsk.pF = pF;
0158 equil_magnetic.eqdsk.pp = pp;
0159 equil_magnetic.eqdsk.pq = pq;
0160 equil_magnetic.eqdsk.Rp = Rp;
0161 equil_magnetic.eqdsk.Zp = Zp;  
0162 equil_magnetic.eqdsk.Ip = Ip;  
0163 equil_magnetic.eqdsk.Rv = Rv;  
0164 equil_magnetic.eqdsk.Zv = Zv;  
0165 equil_magnetic.eqdsk.Bv = Bv;  
0166 equil_magnetic.eqdsk.pFdFdpsi = pFdFdpsi;  
0167 equil_magnetic.eqdsk.pdpdpsi = pdpdpsi;
0168 %
0169 equil.eqdsk = equil_magnetic.eqdsk;
0170 %
0171 %
0172 eval(['save EQUIL_',savename,' equil'])
0173 %
0174 % Verify consistency on pressure between profiles and EFIT
0175 %
0176 if flagprof && display_mode == 2, 
0177     qe = pc_dke_yp;
0178     %
0179     yp = (yTe.*yne + sum(yzTi.*yzni,1))*qe*1e3;%pressure from Tene data
0180     %
0181     figure(13),clf
0182     %
0183     if ~isempty(yrho),
0184         xlab = '\rho';
0185         pdisp = prho;
0186         ydisp = yrho;
0187     elseif ~isempty(yspsin)
0188         xlab = '\psi_n^{1/2}';
0189         pdisp = pspsin;
0190         ydisp = yspsin;
0191     elseif ~isempty(yspsinT)
0192         xlab = '\psi_{nT}^{1/2}';
0193         pdisp = pspsinT;
0194         ydisp = yspsinT;
0195     end  
0196     %
0197     ylab = 'p (hPa)';
0198     xlim = [0,1];
0199     ylim = NaN;
0200     leg = ['T_en_e + T_in_i';
0201            'EFIT           '];
0202     %
0203     graph1D_jd(ydisp,yp/1e5,0,0,'','','',NaN,xlim,ylim,style,marker,color1,width,siz,gca,red,lspace,bspace);
0204     graph1D_jd(pdisp,pp/1e5,0,0,xlab,ylab,tit,leg,xlim,ylim,style,marker,color2,width,siz);
0205     set(gca,'XTick',[0:0.2:1]);
0206     print_jd(p_opt,[savename,'_p_comp']);
0207 end
0208 
0209 
0210