sep4helmex_yp

   Create the separatrix for the toroidal MHD equilibrium solver HELENA (helmex)

   INPUTS

       - xpoint : parameters of the shape of the separatrix
       - ap : minor radius of the plasma at Zp = 0 (m) [1,1] 
       - Rp : major radius of the plasma (m) [1,1] 
       - Zp : vertical shift of the plasma (m) [1,1] 
       - ntheta: number of poloidal points for the separatrix [1,1]
           (default = 65)

   OUTPUTS

       - R_sep0 :  seperatrix major radius [1,ntheta]
       - Z_sep0 :  seperatrix vertical position [1,ntheta]


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

0001 function [R_sep0,Z_sep0] = sep4helmex_yp(xpoint,ntheta)
0002 %
0003 %   Create the separatrix for the toroidal MHD equilibrium solver HELENA (helmex)
0004 %
0005 %   INPUTS
0006 %
0007 %       - xpoint : parameters of the shape of the separatrix
0008 %       - ap : minor radius of the plasma at Zp = 0 (m) [1,1]
0009 %       - Rp : major radius of the plasma (m) [1,1]
0010 %       - Zp : vertical shift of the plasma (m) [1,1]
0011 %       - ntheta: number of poloidal points for the separatrix [1,1]
0012 %           (default = 65)
0013 %
0014 %   OUTPUTS
0015 %
0016 %       - R_sep0 :  seperatrix major radius [1,ntheta]
0017 %       - Z_sep0 :  seperatrix vertical position [1,ntheta]
0018 %
0019 %
0020 %   by Y. Peysson (CEA/DSM/IRFM,yves.peysson@cea.fr) and J. Decker (EPFL/CRPP,joan.decker@epfl.ch)
0021 %
0022 if nargin < 5,
0023     ntheta = 65;
0024 end
0025 %
0026 if nargin < 1,
0027     error('No enough input arguments in sep4helmex_yp.m');
0028 end
0029 %
0030 ap = xpoint(9);
0031 Rp = xpoint(10);
0032 Zp = xpoint(11);
0033 %
0034 lim_top = atan(0.5*xpoint(1)/(1-xpoint(2)))*180/pi;%In order to avoid convergence problem in separatrix MEX file
0035 if xpoint(4) > lim_top,xpoint(4) = lim_top - 0.1;end
0036 lim_bottom = atan(0.5*xpoint(5)/(1-xpoint(6)))*180/pi;%In order to avoid convergence problem in separatrix MEX file
0037 if xpoint(8) > lim_bottom,xpoint(8) = lim_bottom - 0.1;end
0038 %
0039 [x1,y1] = separatrice(xpoint,1);%The LCS is approximated by four ellipses (mode = 1)
0040 %
0041 rr = x1*ap;
0042 zz = y1*ap;
0043 %
0044 alpha  = unwrap(angle(rr+ i .* zz));%poloidal angle
0045 if alpha(end) < -5
0046     alpha = alpha - alpha(end);
0047 end
0048 %
0049 theta_sep  = linspace(0,2*pi,ntheta);% poloidal points
0050 %
0051 ind = find(diff(alpha) == 0);
0052 alpha(ind) = [];
0053 %
0054 rr(ind) = [];
0055 zz(ind) = [];
0056 %
0057 R_sep0 = spline(alpha',rr',theta_sep) + Rp;%resampling
0058 Z_sep0 = spline(alpha',zz',theta_sep) + Zp;%resampling
0059