idealequilcyl_yp

 Create ideal cylindrical magnetic equilibrium with circular concentric flux surfaces from parametric description

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

0001 function [ppsin,ppsi_apRp,theta,x,y,Bx,By,BPHI,pBpp,pq_Rpap,pj,pmag,Ip_test,prho,pspsin,pspsinT,ppsiT] = idealequilcyl_yp(ap,Rp,Zp,Bt,Ip,qmin_Rpap,eq,qopt,npsi,ntheta,lambda)
0002 %
0003 % Create ideal cylindrical magnetic equilibrium with circular concentric flux surfaces from parametric description
0004 %
0005 % by Y.Peysson (CEA/DSM/IRFM, yves.peysson@cea.fr) and J. Decker (CEA/DSM/IRFM, joan.decker@cea.fr)
0006 %
0007 if nargin < 11,
0008     lambda = NaN;% 65 poloidal grid points by default
0009 end
0010 %
0011 if nargin < 10,
0012     ntheta = 65;% 65 poloidal grid points by default
0013 end
0014 %
0015 if nargin < 9,
0016     npsi = 101;% 101 radial profile grid points by default
0017 end
0018 %
0019 if nargin < 8,
0020     qopt = 0;% 101 radial profile grid points by default
0021 end
0022 %
0023 [qe,me,mp,mn,e0,mu0] = pc_dke_yp;
0024 %
0025 Bpa = mu0*Ip*1e6/(2*pi*ap);% Poloidal Field at the edge, Ampere's theorem (circular plasma cross-section). Moreover, in the following expressions, term of order (ap/Rp)^2 or higher are neglected.
0026 qmax_Rpap = abs(Bt/Bpa);%Normalized cylindrical safety factor at the plasma edge
0027 %
0028 prho = linspace(0,1,npsi);%Initial grid
0029 theta = linspace(0,2*pi,ntheta);%Initial grid
0030 %
0031 Xrho = repmat(prho',[1,ntheta]);
0032 Xtheta = repmat(theta,[npsi,1]);
0033 %
0034 x = ap*Xrho.*cos(Xtheta);
0035 y = ap*Xrho.*sin(Xtheta);
0036 %
0037 if isempty(lambda) || isnan(lambda),
0038     if qopt == 2,
0039         if eq >= 1,
0040             pq_Rpap = qmax_Rpap*prho.^2./(1.0 - (1.0 - prho.^2).^(1.0 + eq));%Safety factor profile deduced from a current profile proportional to (1-rho^2)^eq
0041             pq_Rpap(1) = qmax_Rpap/(1.0 + eq);
0042             pdq_Rpap = 2*pq_Rpap.*(1.0 - pq_Rpap*(eq + 1.0).*(1.0 - prho.^2)/qmax_Rpap)./prho;
0043             pdq_Rpap(1) = 0;
0044         else
0045             error('Error: eq >= 1.0 for qopt = 2 !');
0046         end
0047     elseif qopt == 1,
0048          pq_Rpap = (qmax_Rpap - qmin_Rpap)*prho.^eq + qmin_Rpap;%Current profile defined by qmin and eq
0049          pdq_Rpap = eq*(qmax_Rpap - qmin_Rpap)*prho.^(eq-1);
0050     elseif qopt == 0,
0051          pq_Rpap = qmax_Rpap*ones(1,npsi);%Current profile defined with constant q (uniform current density)
0052          pdq_Rpap = zeros(1,npsi);
0053     else
0054         error('This safety factor model is not recognized.')
0055     end
0056     %
0057     pBphi = Bt*ones(size(prho));%
0058     pBpp = sign(Bpa)*prho*abs(Bt)./pq_Rpap;
0059 else
0060     pq_Rpap = prho.*besselj(0,lambda*prho)./besselj(1,lambda*prho);%Safety factor RFP BFM model
0061     pdq_Rpap = besselj(0,lambda*prho)./besselj(1,lambda*prho)...
0062                 - lambda*prho - besselj(0,lambda*prho)./besselj(1,lambda*prho).^2.*(besselj(0,lambda*prho) - lambda*prho.*besselj(2,lambda*prho));
0063     %
0064     pBphi = Bt*besselj(0,lambda*prho);%RFP BFM model
0065     pBpp = Bt*besselj(1,lambda*prho);%RFP BFM model
0066     %
0067     ppsi_apRp_theoric = ap*ap*Bt*(1.0 - besselj(0,lambda*prho))/lambda; 
0068 end
0069 %
0070 %Note: term of order (ap/Rp)^2 or higher are neglected.
0071 %
0072 XBpp = repmat(pBpp',[1,ntheta]);
0073 %
0074 spBpp = triu(repmat((pBpp(1:npsi-1) + pBpp(2:npsi))'/2,[1,npsi-1]));
0075 pdprho = diff(prho);
0076 ppsi_apRp = [0,ap^2*integral_dke_jd(pdprho,spBpp,1)'];%psi_apRp = psi(poloidal usual)*ap/Rp for the cylindrical limit
0077 ppsin = ppsi_apRp/max(ppsi_apRp);%Uniform normalized magnetic flux surface grid
0078 %
0079 pspsin = sqrt(ppsin);
0080 %
0081 ppsiT = pi*Bt*ap^2*prho.^2;%Toroidal magnetic flux
0082 ppsinT = ppsiT/ppsiT(end);%Normalized toroidal magnetic flux
0083 pspsinT = sqrt(ppsinT);
0084 %
0085 Upsilon = 1 + x/Rp;%Toroidal factor. Upsilon = 1 for cylindrical magnetic configuration
0086 %
0087 Bx = -XBpp.*sin(Xtheta)./Upsilon;
0088 By = XBpp.*cos(Xtheta)./Upsilon;
0089 BPHI = repmat(pBphi',[1,ntheta])./Upsilon;
0090 %
0091 pj = (Bt/ap/mu0).*(2./pq_Rpap - pdq_Rpap.*prho./pq_Rpap.^2);%Current density profile (A.m-2)
0092 %
0093 for ii = 1:length(prho),
0094     pIp(ii) = 2*pi*trapz_dke_yp([prho(1:ii)',prho(1:ii)'.*pj(1:ii)'])*ap*ap; 
0095     pmag(ii) = -trapz_dke_yp([prho(1:ii)',pBpp(1:ii)'.*pj(1:ii)'])*ap;
0096 end
0097 pIp(2) = pIp(1);
0098 pmag(2) = pmag(1);
0099 Ip_test = pIp(end)/1e6;%Plasma current (MA)
0100 pmag = (pmag - pmag(end))/qe/1000;%Total pressure consistent with the cylindrical magnetic equilibrium. Pressure is enforced to zero at the plasma edge (special units for helmex77 in keV.m-3)
0101 %