equil_magnetic_ACCOME_yp

 This function reads and returns the magnetic equilibrium parameters from files created by ACCOME or EQDSK for the COMPASS czech tokamak

    INPUTS:

        - filename: equilibrium data and save name
       - jbdirections: current and toroidal B-fields directions (jbdirections.Ip, jbdirections.Bt = 'clockwise' or 'anti-clockwise')
           default: empty
        - display_mode: option for figure display
        - p_opt: option for printing

    OUTPUTS:

        - equil_magnetic: magnetic equilibrium [struct]
       - prho: normalized minor radius [1,nr]
        - psin: normalized poloidal flux [1,nr]

 by Y. Peysson (yves.peysson@cea.fr, CEA/IRFM) and J. Decker (joan.decker@cea.fr, CEA/IRFM)

0001 function [equil_magnetic,prho,psin] = equil_magnetic_ACCOME_yp(filename,jbdirections,display_mode,p_opt,dispname,savename)
0002 %
0003 % This function reads and returns the magnetic equilibrium parameters from files created by ACCOME or EQDSK for the COMPASS czech tokamak
0004 %
0005 %    INPUTS:
0006 %
0007 %        - filename: equilibrium data and save name
0008 %       - jbdirections: current and toroidal B-fields directions (jbdirections.Ip, jbdirections.Bt = 'clockwise' or 'anti-clockwise')
0009 %           default: empty
0010 %        - display_mode: option for figure display
0011 %        - p_opt: option for printing
0012 %
0013 %    OUTPUTS:
0014 %
0015 %        - equil_magnetic: magnetic equilibrium [struct]
0016 %       - prho: normalized minor radius [1,nr]
0017 %        - psin: normalized poloidal flux [1,nr]
0018 %
0019 % by Y. Peysson (yves.peysson@cea.fr, CEA/IRFM) and J. Decker (joan.decker@cea.fr, CEA/IRFM)
0020 %
0021 if nargin < 6,
0022     dispname = filename;
0023 end
0024 %
0025 if nargin < 5,
0026     dispname = filename;
0027 end
0028 %
0029 if nargin < 4,
0030     p_opt = -1;
0031 end
0032 %
0033 if nargin < 3,
0034     display_mode = 0;
0035 end
0036 %
0037 if nargin < 2,
0038     jbdirections = [];
0039 end
0040 %
0041 method = 'spline';%'linear';
0042 nt = 65;
0043 %
0044 % Reads ACCOME data
0045 %
0046 if exist([filename,'.bfield'],'file'),
0047     disp(' ---> ACCOME toroidal MHD equilibrium file format found.');
0048     [x,y,xypsin,psia,rpsin,xyBPHI,xyBR,xyBZ,Rp,Zp] = accome_yp([filename,'.bfield'],display_mode,p_opt,dispname,savename);
0049 elseif exist([filename,'.eqdsk'],'file'),
0050     disp(' ---> EQDSK toroidal MHD equilibrium file format found.');
0051     [x,y,xypsin,psia,rpsin,rF,rp,rq,Rp,Zp,Ip,Rv,Zv,Bv,rFdFdpsi,rdpdpsi] = eqdsk_jd([filename,'.eqdsk'],display_mode,p_opt,dispname,savename);
0052 else
0053     error(' ---> Unknown toroidal MHD equilibrium file format found.')
0054 end
0055 %
0056 if ~isempty(jbdirections),
0057     %
0058     if strcmp(jbdirections.Ip,'clockwise')
0059         if psia > 0,
0060             if exist('xyBR','var') && exist('xyBZ','var')
0061                 if xyBR(end,end-1) < 0 && xyBZ(end,end-1) > 0
0062                     disp(' ---> The input magnetic equilibrium is consistent with Ip clockwise: LUKE signs convention)')
0063                 else
0064                     error('Signs of psia, xyBR, xyBZ in the input magnetic equilibrium are inconsistent')
0065                 end
0066             end 
0067         else
0068            psia = -psia; 
0069            disp(' ---> WARNING: psia positive enforced, Ip positive enforced (clockwise: LUKE signs convention)') 
0070             if exist('xyBR','var') && exist('xyBZ','var')
0071                 if xyBR(end,end-1) > 0 && xyBZ(end,end-1) < 0
0072                     xyBR = - xyBR;
0073                     xyBZ = - xyBZ;
0074                 else
0075                     error('Signs of psia, xyBR, xyBZ in the input magnetic equilibrium are inconsistent')
0076                 end
0077             end
0078         end
0079     else
0080         if psia < 0,
0081             if exist('xyBR','var') && exist('xyBZ','var')
0082                 if xyBR(end,end-1) > 0 && xyBZ(end,end-1) < 0
0083                     disp(' ---> The input magnetic equilibrium is consistent with Ip anti-clockwise: LUKE signs convention)')
0084                 else
0085                     error('Signs of psia, xyBR, xyBZ in the input magnetic equilibrium are inconsistent')
0086                 end
0087             end 
0088         else
0089            psia = -psia; 
0090            disp(' ---> WARNING: psia negtaive enforced, Ip negative enforced (counter-clockwise: LUKE signs convention)') 
0091             if exist('xyBR','var') && exist('xyBZ','var')
0092                 if xyBR(end,end-1) < 0 && xyBZ(end,end-1) > 0
0093                     xyBR = - xyBR;
0094                     xyBZ = - xyBZ;
0095                 else
0096                     error('Signs of psia, xyBR, xyBZ in the input magnetic equilibrium are inconsistent')
0097                 end
0098             end
0099         end     
0100     end
0101     %
0102     if exist('Ip','var'),Ip = abs(Ip)*sign(psia);end
0103     %
0104     if strcmp(jbdirections.Bt,'clockwise'),
0105         disp(' ---> WARNING: BPHI positive enforced (clockwise: LUKE signs convention)')
0106         if exist('xyBPHI','var'), xyBPHI = abs(xyBPHI);end
0107         if exist('rF','var'),rF = abs(rF);end 
0108     else
0109         disp(' ---> WARNING: BPHI negative enforced (counter-clockwise: LUKE signs convention)')
0110         if exist('xyBPHI','var'), xyBPHI = -abs(xyBPHI);end
0111         if exist('rF','var'),rF = -abs(rF);end
0112     end
0113     %
0114 else
0115     disp(' ---> WARNING: LUKE signs conventions not enforced for Ip and Bt, check carefully the data consistency in ''equil'' structure.')
0116 end
0117 %
0118 % Transforms on a (psi,theta) grid
0119 %
0120 [ptx,pty,ptdpsindx,ptdpsindy,psin,theta] = rz2pt_jd(x,y,xypsin,rpsin,nt,method,display_mode,p_opt,savename);
0121 %
0122 if exist('xyBPHI','var');
0123     ptBPHI = griddata(x,y,xyBPHI',ptx,pty);
0124 end
0125 %
0126 np = length(psin);
0127 ap = ptx(np,1);%Plasma minor radius on the LFS midplane
0128 psi_apRp = psin*psia*ap/Rp;%The poloidal flux coordinate is normalized by the aspect ratio
0129 ptR = Rp + ptx;
0130 %
0131 % Note that in DKE we have B = F.gradphi + gradphi x gradpsi.
0132 %
0133 %WARNING: not very accurate because psi grid provided by ACCOME seems to be noisy in the vicinity of the plasma center
0134 %
0135 if ~exist('xyBR','var') || ~exist('xyBZ','var') || isempty(xyBR) || isempty(xyBZ),
0136     %
0137     ptBx = -ptdpsindy*psia./ptR;
0138     ptBy = +ptdpsindx*psia./ptR;
0139     %
0140     ptF = repmat(rF.',[1,nt]);
0141     ptBPHI = ptF./ptR;
0142     %
0143     if exist([filename,'.bfield'],'file'),
0144         disp('***********************************************************************************************************');
0145         disp('WARNING: The BR and BZ magnetic field components are not given explicitely by ACCOME (not in the source file)');
0146         disp('         BR and BZ are calculated from the derivatives of the poloidal magnetic flux.');
0147         disp('         BR and BZ are likely not very accurate in the vicinity of the magnetic axis.');
0148         disp('         It is recommanded to re-create if possible an ACCOME source file for LUKE.');
0149         disp('***********************************************************************************************************');
0150     end
0151 else
0152     ptBx = griddata(x,y,xyBR',ptx,pty);
0153     ptBy = griddata(x,y,xyBZ',ptx,pty);
0154 end
0155 %
0156 ptBP = sqrt(ptBx.^2 +ptBy.^2);
0157 ptB = sqrt(ptBP.^2 + ptBPHI.^2);
0158 %
0159 prho = ptx(:,1).'/ap;
0160 %
0161 if display_mode == 2, 
0162 %
0163     rx = [-fliplr(prho(2:np)),prho];
0164     rBx = [fliplr(ptBx(2:np,(nt+1)/2)'),ptBx(:,1)'];
0165     rBy = [fliplr(ptBy(2:np,(nt+1)/2)'),ptBy(:,1)'];
0166     rBPHI = [fliplr(ptBPHI(2:np,(nt+1)/2)'),ptBPHI(:,1)'];
0167     rBP = sqrt(rBx.^2 + rBy.^2);
0168     rB = sqrt(rBP.^2 + rBPHI.^2);
0169 %
0170     style = '-';
0171     marker = 'none';
0172     color1 = 'r';
0173     color2 = 'b';
0174     color3 = 'k';
0175     width = 2;
0176     siz = 20;
0177     red = 0.9;
0178     lspace = 0.7;
0179     bspace = 0.6;
0180     tit = [dispname];
0181     %
0182     figure('Name','Poloidal and toroidal magnetic fields'),clf
0183     %
0184     xlab = 'r/a';
0185     ylab = 'B(T)';
0186     xlim = [-1,1];
0187     ylim = NaN;
0188     leg = ['B_T';'B_P';'B  '];
0189 %
0190     graph1D_jd(rx,abs(rBPHI),0,0,'','','',NaN,xlim,ylim,style,marker,color1,width,siz,gca,red,lspace,bspace);
0191     graph1D_jd(rx,rBP,0,0,'','','',NaN,xlim,ylim,style,marker,color2,width);
0192     graph1D_jd(rx,rB,0,0,xlab,ylab,tit,leg,xlim,ylim,style,marker,color3,width);
0193     print_jd(p_opt,[savename,'_Bfield']);
0194 %
0195     figure('Name','Total magnetic fields'),clf
0196     graph1D_jd(theta,ptB(2,:),0,0,'','','',NaN,'','',style,marker,'b',width);hold on
0197     plot(theta,ptB([2:end],:))
0198     legend(['r/a = ',num2str(prho(2))])
0199     axis([0,2*pi,0,ceil(max(max(ptB([2:end],:)))*10)/10])
0200     xlabel('\theta (rd)'),ylabel('B(T)'),title([dispname])
0201 end
0202 %
0203 ptx(:,65) = ptx(:,1);
0204 pty(:,65) = pty(:,1);
0205 ptBx(:,65) = ptBx(:,1);
0206 ptBy(:,65) = ptBy(:,1);
0207 ptBPHI(:,65) = ptBPHI(:,1);
0208 %
0209 % equil structure
0210 %
0211 equil_magnetic.psi_apRp = psi_apRp;
0212 equil_magnetic.theta = theta;
0213 equil_magnetic.ptx = ptx;
0214 equil_magnetic.pty = pty;
0215 equil_magnetic.ptBx = ptBx;
0216 equil_magnetic.ptBy = ptBy;
0217 equil_magnetic.ptBPHI = ptBPHI;
0218 equil_magnetic.Rp = Rp;
0219 equil_magnetic.Zp = Zp;
0220 %
0221 if exist([filename,'.eqdsk'],'file'),
0222     equil_magnetic.eqdsk.x = x;
0223     equil_magnetic.eqdsk.y = y;
0224     equil_magnetic.eqdsk.xypsin = xypsin;
0225     equil_magnetic.eqdsk.psia = psia;
0226     equil_magnetic.eqdsk.psin = psin;
0227     equil_magnetic.eqdsk.pF = rF;
0228     equil_magnetic.eqdsk.pp = rp;
0229     equil_magnetic.eqdsk.pq = rq;
0230     equil_magnetic.eqdsk.Rp = Rp;
0231     equil_magnetic.eqdsk.Zp = Zp;  
0232     equil_magnetic.eqdsk.Ip = Ip;  
0233     equil_magnetic.eqdsk.Rv = Rv;  
0234     equil_magnetic.eqdsk.Zv = Zv;  
0235     equil_magnetic.eqdsk.Bv = Bv;  
0236     equil_magnetic.eqdsk.pFdFdpsi = rFdFdpsi;  
0237     equil_magnetic.eqdsk.pdpdpsi = rdpdpsi;  
0238 end
0239 %