equil_magnetic_EFIT_jd

 This function reads and returns the magnetic equilibrium parameters from 
 eqdsk files created
 by EFIT and saved in the formatted text file "filename.efit".

    INPUTS:

        - eqdsk: equilibrium data and save name (if savename not
        specified). If eqdsk is a structure, no need of file in eqdsk_jd
        - display_mode: option for figure display
        - p_opt: option for printing

    OUTPUTS:

       - equil_magnetic: LUKE magnetic equilibrium structure
       - prho: normalized radial coordinate (geometrical as defined in LUKE at theta_pol = 0)
       - pspsin: normalized radial coordinate, sqrt(normalized poloidal flux)
       - pspsinT: normalized radial coordinate,  sqrt(normalized toroidal flux)
       - pp: pressure values at nnr equispaced points in psi from psimag to psilim (see eqdsk_jd.m)

 by J. Decker (joan.decker@cea.fr) 11/01/04

0001 function [equil_magnetic,prho,pspsin,pspsinT,pp] = equil_magnetic_EFIT_jd(eqdsk,display_mode,p_opt,dispname,savename,opt)
0002 %
0003 % This function reads and returns the magnetic equilibrium parameters from
0004 % eqdsk files created
0005 % by EFIT and saved in the formatted text file "filename.efit".
0006 %
0007 %    INPUTS:
0008 %
0009 %        - eqdsk: equilibrium data and save name (if savename not
0010 %        specified). If eqdsk is a structure, no need of file in eqdsk_jd
0011 %        - display_mode: option for figure display
0012 %        - p_opt: option for printing
0013 %
0014 %    OUTPUTS:
0015 %
0016 %       - equil_magnetic: LUKE magnetic equilibrium structure
0017 %       - prho: normalized radial coordinate (geometrical as defined in LUKE at theta_pol = 0)
0018 %       - pspsin: normalized radial coordinate, sqrt(normalized poloidal flux)
0019 %       - pspsinT: normalized radial coordinate,  sqrt(normalized toroidal flux)
0020 %       - pp: pressure values at nnr equispaced points in psi from psimag to psilim (see eqdsk_jd.m)
0021 %
0022 % by J. Decker (joan.decker@cea.fr) 11/01/04
0023 %
0024 if nargin < 6,
0025     opt.sfac = 10;% factor for (x,y) supergrid
0026     opt.axenforce = true;% option to enforce center position and psi value in supergrid
0027 end
0028 %
0029 if nargin < 5,
0030     savename = 'equil';
0031 end
0032 %
0033 if nargin < 4,
0034     dispname = '';
0035 end
0036 %
0037 if nargin < 3,
0038     p_opt = -1;
0039 end
0040 %
0041 if nargin < 2,
0042     display_mode = 0;
0043 end
0044 %
0045 method = 'spline';%'linear';
0046 nt = 65;
0047 %
0048 % Reads EFIT eqdsk data
0049 %
0050 [x,y,xypsin,psia,psin,pF,pp,pq,Rp,Zp,Ip,Rv,Zv,Bv,pFdFdpsi,pdpdpsi] = eqdsk_jd(eqdsk,display_mode,p_opt,dispname,savename);
0051 %
0052 % Transforms on a (psi,theta) grid
0053 %
0054 [ptx,pty,ptdpsindx,ptdpsindy,psin,theta] = rz2pt_jd(x,y,xypsin,psin,nt,method,display_mode,p_opt,savename,opt);
0055 %
0056 np = length(psin);
0057 ap = ptx(np,1);%Plasma minor radius on the LFS midplane
0058 psi_apRp = psin*psia*ap/Rp;%The poloidal flux coordinate is normalized by the aspect ratio
0059 ptR = Rp + ptx;
0060 %
0061 % Note that in DKE we have B = F.gradphi + gradphi x gradpsi.
0062 %
0063 ptBx = -ptdpsindy*psia./ptR;
0064 ptBy = ptdpsindx*psia./ptR;
0065 %
0066 ptF = repmat(pF.',[1,nt]);
0067 ptBPHI = ptF./ptR;
0068 %
0069 ptBP = sqrt(ptBx.^2 +ptBy.^2);
0070 ptB = sqrt(ptBP.^2 + ptBPHI.^2);
0071 %
0072 prho = ptx(:,1).'/ap;
0073 %
0074 pq_h = (pq(1:end-1) + pq(2:end))/2;
0075 pdpsin = diff(psin);
0076 ppsiT = [0,cumsum(pq_h.*pdpsin)];
0077 pspsinT = sqrt(ppsiT/ppsiT(end));
0078 %
0079 pspsin = sqrt(psin);
0080 %
0081 if display_mode == 2, 
0082 %
0083     rx = [-fliplr(prho(2:np)),prho];
0084     rBx = [fliplr(ptBx(2:np,(nt+1)/2)'),ptBx(:,1)'];
0085     rBy = [fliplr(ptBy(2:np,(nt+1)/2)'),ptBy(:,1)'];
0086     rBPHI = [fliplr(ptBPHI(2:np,(nt+1)/2)'),ptBPHI(:,1)'];
0087     rBP = sqrt(rBx.^2 + rBy.^2);
0088     rB = sqrt(rBP.^2 + rBPHI.^2);
0089 %
0090     style = '-';
0091     marker = 'none';
0092     color1 = 'r';
0093     color2 = 'b';
0094     color3 = 'k';
0095     width = 2;
0096     siz = 20;
0097     red = 0.9;
0098     lspace = 0.7;
0099     bspace = 0.6;
0100     tit = [dispname];
0101     %
0102     figure(12),clf
0103     %
0104     xlab = 'r/a';
0105     ylab = 'B(T)';
0106     xlim = [-1,1];
0107     ylim = NaN;
0108     leg = ['B_T';'B_P';'B  '];
0109 %
0110     graph1D_jd(rx,abs(rBPHI),0,0,'','','',NaN,xlim,ylim,style,marker,color1,width,siz,gca,red,lspace,bspace);
0111     graph1D_jd(rx,rBP,0,0,'','','',NaN,xlim,ylim,style,marker,color2,width);
0112     graph1D_jd(rx,rB,0,0,xlab,ylab,tit,leg,xlim,ylim,style,marker,color3,width);
0113     print_jd(p_opt,[savename,'_Bfield']);
0114 %
0115     figure(13),clf
0116     graph1D_jd(theta,ptB(2,:),0,0,'','','',NaN,'','',style,marker,'b',width);hold on
0117     plot(theta,ptB([2:end],:))
0118     legend(['r/a = ',num2str(prho(2))])
0119     axis([0,2*pi,0,ceil(max(max(ptB([2:end],:)))*10)/10])
0120     xlabel('\theta (rd)'),ylabel('B(T)'),title([dispname])
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_magnetic.psi_apRp = psi_apRp;
0132 equil_magnetic.theta = theta;
0133 equil_magnetic.ptx = ptx;
0134 equil_magnetic.pty = pty;
0135 equil_magnetic.ptBx = ptBx;
0136 equil_magnetic.ptBy = ptBy;
0137 equil_magnetic.ptBPHI = ptBPHI;
0138 equil_magnetic.Rp = Rp;
0139 equil_magnetic.Zp = Zp;
0140 equil_magnetic.equi_src = eqdsk;
0141 %
0142 equil_magnetic.eqdsk.x = x;
0143 equil_magnetic.eqdsk.y = y;
0144 equil_magnetic.eqdsk.xypsin = xypsin;
0145 equil_magnetic.eqdsk.psia = psia;
0146 equil_magnetic.eqdsk.psin = psin;
0147 equil_magnetic.eqdsk.pF = pF;
0148 equil_magnetic.eqdsk.pp = pp;
0149 equil_magnetic.eqdsk.pq = pq;
0150 equil_magnetic.eqdsk.Rp = Rp;
0151 equil_magnetic.eqdsk.Zp = Zp;  
0152 equil_magnetic.eqdsk.Ip = Ip;  
0153 equil_magnetic.eqdsk.Rv = Rv;  
0154 equil_magnetic.eqdsk.Zv = Zv;  
0155 equil_magnetic.eqdsk.Bv = Bv;  
0156 equil_magnetic.eqdsk.pFdFdpsi = pFdFdpsi;  
0157 equil_magnetic.eqdsk.pdpdpsi = pdpdpsi;  
0158 %