equildipole_fxd

 Create dipole magnetic equilibrium

 by F.-X Duthoit (CEA/DSM/IRFM, francois-xavier.duthoit@cea.fr),
    J. Decker (CEA/DSM/IRFM, joan.decker@cea.fr) and
    Y. Peysson (CEA/DSM/IRFM, yves.paysson@cea.fr)

0001 function equil_magnetic = equildipole_fxd(ap,RT,M,prho,ntheta,test_mode,method)
0002 %
0003 % Create dipole magnetic equilibrium
0004 %
0005 % by F.-X Duthoit (CEA/DSM/IRFM, francois-xavier.duthoit@cea.fr),
0006 %    J. Decker (CEA/DSM/IRFM, joan.decker@cea.fr) and
0007 %    Y. Peysson (CEA/DSM/IRFM, yves.paysson@cea.fr)
0008 %
0009 equil_magnetic = [];
0010 return
0011 %TODO : fix error that block LUKE install YP 20/01/17
0012 %
0013 if nargin < 7,
0014     method = 1;
0015 end
0016 %
0017 if nargin < 6,
0018     test_mode = 0;
0019 end
0020 %
0021 if nargin < 5,
0022     ntheta = 65;
0023 end
0024 %
0025 if nargin < 4,
0026     prho = linspace(0,1,101); 
0027 end
0028 %
0029 if (ntheta - 1)/4 ~= round((ntheta - 1)/4),
0030     error('ntheta - 1 must be a multiple of 4')
0031 end
0032 %
0033 it1 = (ntheta - 1)/4 + 1;
0034 it2 = (ntheta - 1)/2 + 1;
0035 it3 = 3*(ntheta - 1)/4 + 1;
0036 %
0037 iar = ap/RT;
0038 %
0039 [qe,me,mp,mn,e0,mu0] = pc_dke_yp;
0040 %
0041 npsi = length(prho);%radial grid rho = (R - RT)/ap taken on LFS horizontal midplane
0042 theta = linspace(0,2*pi,ntheta);%poloidal grid
0043 %
0044 pR0 = 1 + iar*prho;
0045 %
0046 if method == 0,
0047     %
0048     nlambda = 1001;
0049     N = fix((ntheta - 1)/2);
0050     %
0051     lambda = linspace(0,pi,nlambda);%vertical angle in spherical view
0052     %
0053     XR0 = repmat(pR0,[nlambda,1]);
0054     Xlambda = repmat(lambda.',[1,npsi]);
0055     %
0056     Xrs = XR0.*sin(Xlambda).^2;
0057     %
0058     Xx = Xrs.*sin(Xlambda) - 1;
0059     Xy = Xrs.*cos(Xlambda);
0060     %
0061     Xtheta = atan(Xy./Xx) + pi*(Xx < 0).*(Xy >= 0) - pi*(Xx < 0).*(Xy < 0);% -pi < theta < pi
0062     Ydtheta = -diff(Xtheta);
0063     Ytheta = (Xtheta(1:end-1,:) + Xtheta(2:end,:))/2;
0064     Yx = (Xx(1:end-1,:) + Xx(2:end,:))/2;
0065     Yy = (Xy(1:end-1,:) + Xy(2:end,:))/2;
0066     %
0067     n = 1:N;
0068     %
0069     Ynx = repmat(Yx,[1,1,N]);
0070     Yny = repmat(Yy,[1,1,N]);
0071     Yntheta = repmat(Ytheta,[1,1,N]);
0072     Yndtheta = repmat(Ydtheta,[1,1,N]);
0073     %
0074     Ynn = repmat(reshape(n,[1,1,N]),[nlambda-1,npsi,1]);
0075     %
0076     % Note that x is even is theta while y is odd
0077     %
0078     npan = reshape(sum(Ynx.*cos(Ynn.*Yntheta).*Yndtheta)/pi,[npsi,N]).';
0079     pa0 = sum(Yx.*Ydtheta)/pi;
0080     npbn = reshape(sum(Yny.*sin(Ynn.*Yntheta).*Yndtheta)/pi,[npsi,N]).';
0081     %
0082     nptan = repmat(npan,[1,1,ntheta]);
0083     nptbn = repmat(npbn,[1,1,ntheta]);
0084     nptn = repmat(n.',[1,npsi,ntheta]);
0085     nptheta = repmat(reshape(theta,[1,1,ntheta]),[N,npsi,1]);
0086     pta0 = repmat(pa0.',[1,ntheta]);
0087     %
0088     ptx = reshape(sum(nptan.*cos(nptn.*nptheta)),[npsi,ntheta]) + pta0/2;
0089     pty = reshape(sum(nptbn.*sin(nptn.*nptheta)),[npsi,ntheta]);
0090     %
0091     ptR = ptx + 1;
0092     ptZ = pty;
0093     %
0094     ptrhos = sqrt(ptR.^2+ptZ.^2);
0095     %
0096 else
0097     %
0098     ptR = NaN(npsi,ntheta);
0099     %
0100     for it = 1:ntheta,
0101         %
0102         c = cos(theta(it));
0103         c6 = c^6;
0104         s2 = sin(theta(it))^2;
0105         s4 = s2*s2;
0106         s6 = s4*s2;
0107         %
0108         A3 = 3*s2*(1+4*s2);
0109         A = [1,-6*s2,NaN,-4*s4*(3+2*s2),3*s4*(1+4*s2),-6*s6,s6];
0110         %
0111         for ir = 2:npsi,
0112             %
0113             A(3) = A3 - c6*pR0(ir)*pR0(ir);
0114             %
0115             R = roots(A);
0116             %
0117             iR = find(imag(R) == 0 & real(R) >= 0 & (real(R) - 1)*c >= 0);
0118             %
0119             if isempty(iR),
0120                 error('no root found')
0121             elseif length(iR) > 1 && max(abs(diff(R(iR)))) > 0,
0122                 error('no unique root found')
0123             end
0124             %
0125             ptR(ir,it) = R(iR(1));
0126             %
0127         end
0128     end
0129     %
0130     % particular case rho = 0
0131     %
0132     for it = 1:ntheta,
0133         %
0134         c = cos(theta(it));
0135         %
0136         if c < 0,
0137             c6 = c^6;
0138             s2 = sin(theta(it))^2;
0139             s4 = s2*s2;
0140             s6 = s4*s2;
0141             %
0142             A = [1,-6*s2,3*s2*(1+4*s2)-c6,-4*s4*(3+2*s2),3*s4*(1+4*s2),-6*s6,s6];
0143             %
0144             R = roots(A);
0145             %
0146             iR = find(imag(R) == 0 & real(R) >= 0);
0147             %
0148             if isempty(iR),
0149                 error('no root found')
0150             end
0151             %
0152             ptR(1,it) = min(R(iR));
0153             %
0154         else
0155             %
0156             ptR(1,it) = 1;
0157             %
0158         end
0159     end
0160     %
0161     ptR(:,1) = pR0;% to avoid roundup errors
0162     ptR(:,it1) = 1;% to avoid roundup errors
0163     ptR(:,it2) = 0;% to avoid roundup errors
0164     ptR(:,it3) = 1;% to avoid roundup errors
0165     ptR(:,end) = pR0;% to avoid roundup errors
0166     %
0167     ptR0 = repmat(pR0.',[1,ntheta]);
0168     ptrhos = (ptR0.*ptR.*ptR).^(1/3);
0169     ptZ = sqrt(ptrhos.^2 - ptR.^2).*repmat(sign(sin(theta)),[npsi,1]);
0170     %
0171     ptZ(:,1) = 0;% to avoid roundup errors
0172     ptZ(:,end) = 0;% to avoid roundup errors
0173     %
0174     ptx = ptR - 1;
0175     pty = ptZ;
0176     %
0177 end
0178 %
0179 C = mu0*M/(4*pi);
0180 %
0181 ptBx = 3*C*ptR.*ptZ./ptrhos.^5;
0182 ptBy = C*(2*ptZ.^2-ptR.^2)./ptrhos.^5;
0183 ptBPHI = zeros(npsi,ntheta);
0184 %
0185 ptBx(:,it2) = 0;
0186 ptBy(:,it2) = -sign(M)*Inf;
0187 %
0188 psi = C*(1./pR0 - 1);
0189 psi_apRp = psi*iar;
0190 %
0191 if test_mode
0192     %
0193     % The precision of interpolation is tested by ensuring that psi is constant along the field lines
0194     %
0195     figure(1),
0196     plot(ptx.',pty.');axis('equal');
0197     %
0198     ptmask = ptrhos >= 1;
0199     pntheta = sum(ptmask.');
0200     %
0201     ptpsi = C*(ptR.^2./ptrhos.^3 - 1);
0202     ptpsi(:,it2) = 0;%theta = pi
0203     %
0204     pavg_ptpsi = sum(ptpsi.'.*ptmask.')./pntheta;
0205     pdev_ptpsi = sqrt(sum((ptpsi.' - repmat(pavg_ptpsi,[ntheta,1])).^2.*ptmask.'))./pntheta;
0206     prec_ptpsi = max(abs(ptpsi.' - repmat(psi,[ntheta,1])).*ptmask.');
0207     %
0208     figure(2),
0209     semilogy(prho,abs(pavg_ptpsi - psi)./psi,prho,pdev_ptpsi./psi,prho,prec_ptpsi./psi);
0210     %
0211     phtBV = (ptBy(1:end-1,1) + ptBy(2:end,1))/2;
0212     phtdA_2pi = diff(ptx(:,1)).*(ptR(1:end-1,1) + ptR(2:end,1))/2;
0213     psiI = [0,cumsum(phtBV.*phtdA_2pi).'];
0214     %
0215     figure(3),clf
0216     plot(prho,psi,prho,psiI,'r--')
0217     %
0218 end
0219 %
0220 % renormalization
0221 %
0222 equil_magnetic.psi_apRp = psi_apRp/RT;
0223 equil_magnetic.theta = theta;
0224 equil_magnetic.ptx = ptx*RT;
0225 equil_magnetic.pty = pty*RT;
0226 equil_magnetic.ptBx = ptBx/RT^3;
0227 equil_magnetic.ptBy = ptBy/RT^3;
0228 equil_magnetic.ptBPHI = ptBPHI/RT^3;
0229 %
0230 equil_magnetic.Rp = RT;
0231 equil_magnetic.Zp = 0;
0232 %