equil_magnetic_rz2pt_jd

 This function reads and returns the magnetic equilibrium parameters from 
 magnetic field data on a RZ grid

    INPUTS:

        - equil_magnetic_data: equilibrium data
        - display_mode: option for figure display
        - p_opt: option for printing

 by J. Decker (joan.decker@cea.fr) and Y. Peysson

0001 function [equil_magnetic,prho,pspsin,pspsinT] = equil_magnetic_rz2pt_jd(equil_magnetic_data,display_mode,p_opt,dispname,savename,np,nt,opt_in)
0002 %
0003 % This function reads and returns the magnetic equilibrium parameters from
0004 % magnetic field data on a RZ grid
0005 %
0006 %    INPUTS:
0007 %
0008 %        - equil_magnetic_data: equilibrium data
0009 %        - display_mode: option for figure display
0010 %        - p_opt: option for printing
0011 %
0012 % by J. Decker (joan.decker@cea.fr) and Y. Peysson
0013 %
0014 if nargin < 7,
0015     nt = 65;
0016 end
0017 %
0018 if nargin < 6,
0019     np = 101;
0020 end
0021 %
0022 if nargin < 5,
0023     savename = '';
0024 end
0025 %
0026 if nargin < 4,
0027     dispname = '';
0028 end
0029 %
0030 if nargin < 3,
0031     p_opt = -1;
0032 end
0033 %
0034 if nargin < 2,
0035     display_mode = 0;
0036 end
0037 %
0038 method = 'spline';%'linear';
0039 %
0040 x = equil_magnetic_data.x;
0041 y = equil_magnetic_data.y;
0042 Rp = equil_magnetic_data.Rp;
0043 Zp = equil_magnetic_data.Zp;
0044 xypsin = equil_magnetic_data.xypsin;
0045 psia = equil_magnetic_data.psia;
0046 xBt = equil_magnetic_data.xBt;
0047 ypsin = equil_magnetic_data.psin;
0048 yq = equil_magnetic_data.pq;
0049 %
0050 psin = linspace(0,1,np);
0051 %
0052 if ~struct(opt_in),
0053     opt.sfac = opt_in;
0054     opt.axenforce = true;% option to enforce center position and psi value in supergrid
0055 end
0056 %
0057 % Transforms on a (psi,theta) grid
0058 %
0059 [ptx,pty,ptdpsindx,ptdpsindy,psin,theta] = rz2pt_jd(x,y,xypsin,psin,nt,method,display_mode,p_opt,savename,opt);
0060 %
0061 np = length(psin);
0062 ap = ptx(np,1);%Plasma minor radius on the LFS midplane
0063 psi_apRp = psin*psia*ap/Rp;%The poloidal flux coordinate is normalized by the aspect ratio
0064 ptR = Rp + ptx;
0065 %
0066 % Note that in DKE we have B = F.gradphi + gradphi x gradpsi.
0067 %
0068 ptBx = -ptdpsindy*psia./ptR;
0069 ptBy = ptdpsindx*psia./ptR;
0070 %
0071 ptBPHI = interp1(x,xBt,ptx,method);
0072 %
0073 ptBP = sqrt(ptBx.^2 +ptBy.^2);
0074 ptB = sqrt(ptBP.^2 + ptBPHI.^2);
0075 %
0076 prho = ptx(:,1).'/ap;
0077 %
0078 pq = interp1(ypsin,yq,psin,method,'extrap');
0079 %
0080 pq_h = (pq(1:end-1) + pq(2:end))/2;
0081 pdpsin = diff(psin);
0082 ppsiT = [0,cumsum(pq_h.*pdpsin)];
0083 pspsinT = sqrt(ppsiT/ppsiT(end));
0084 %
0085 pspsin = sqrt(psin);
0086 %
0087 if display_mode == 2, 
0088 %
0089     rx = [-fliplr(prho(2:np)),prho];
0090     rBx = [fliplr(ptBx(2:np,(nt+1)/2)'),ptBx(:,1)'];
0091     rBy = [fliplr(ptBy(2:np,(nt+1)/2)'),ptBy(:,1)'];
0092     rBPHI = [fliplr(ptBPHI(2:np,(nt+1)/2)'),ptBPHI(:,1)'];
0093     rBP = sqrt(rBx.^2 + rBy.^2);
0094     rB = sqrt(rBP.^2 + rBPHI.^2);
0095 %
0096     style = '-';
0097     marker = 'none';
0098     color1 = 'r';
0099     color2 = 'b';
0100     color3 = 'k';
0101     width = 2;
0102     siz = 20;
0103     red = 0.9;
0104     lspace = 0.7;
0105     bspace = 0.6;
0106     tit = [dispname];
0107     %
0108     figure(12),clf
0109     %
0110     xlab = 'r/a';
0111     ylab = 'B(T)';
0112     xlim = [-1,1];
0113     ylim = NaN;
0114     leg = ['B_T';'B_P';'B  '];
0115 %
0116     graph1D_jd(rx,abs(rBPHI),0,0,'','','',NaN,xlim,ylim,style,marker,color1,width,siz,gca,red,lspace,bspace);
0117     graph1D_jd(rx,rBP,0,0,'','','',NaN,xlim,ylim,style,marker,color2,width);
0118     graph1D_jd(rx,rB,0,0,xlab,ylab,tit,leg,xlim,ylim,style,marker,color3,width);
0119     print_jd(p_opt,[savename,'_Bfield']);
0120 %
0121     figure(13),clf
0122     graph1D_jd(theta,ptB(2,:),0,0,'','','',NaN,'','',style,marker,'b',width);hold on
0123     plot(theta,ptB([2:end],:))
0124     legend(['r/a = ',num2str(prho(2))])
0125     axis([0,2*pi,0,ceil(max(max(ptB([2:end],:)))*10)/10])
0126     xlabel('\theta (rd)'),ylabel('B(T)'),title([dispname])
0127 end
0128 %
0129 ptx(:,65) = ptx(:,1);
0130 pty(:,65) = pty(:,1);
0131 ptBx(:,65) = ptBx(:,1);
0132 ptBy(:,65) = ptBy(:,1);
0133 ptBPHI(:,65) = ptBPHI(:,1);
0134 %
0135 % equil structure
0136 %
0137 equil_magnetic.psi_apRp = psi_apRp;
0138 equil_magnetic.theta = theta;
0139 equil_magnetic.ptx = ptx;
0140 equil_magnetic.pty = pty;
0141 equil_magnetic.ptBx = ptBx;
0142 equil_magnetic.ptBy = ptBy;
0143 equil_magnetic.ptBPHI = ptBPHI;
0144 equil_magnetic.Rp = Rp;
0145 equil_magnetic.Zp = Zp;
0146 %