fitequil_yp

 C3PO - Build vectorial tokamak magnetic equilibrium from its numerical counterpart

 Build vectorial tokamak magnetic equilibrium from its numerical counterpart

 Two operating modes :
 1) Spline-Fourier description for (psi,theta) magnetic equilibrium (HELENA for example)

    1-D: use interp1.m MatLab built-in function (Hermite polynomials (pchip), spline, linear,  nearest,...)
    2-D: use Fourier series expansion for the poloidal dependence and the
    interp1.m for the Fourier series coefficients (Hermite polynomials (pchip), spline, linear,
    nearest,...)

 2) Spline-Spline description for (x,y) magnetic equilibrium (EFIT for example, or for TORPEX case)

    1-D: use interp1.m MatLab built-in function (Hermite polynomials (pchip), spline, linear, nearest,...)
    2-D: use interp1.m for the x-direction and the interp1.m again for the coefficients of the 
    previous spline (Hermite polynomials (pchip), spline, bi-cubic, linear, nearest,...)

 By Yves Peysson (CEA-DRFC, yves.peysson@cea.fr) and Joan Decker (CEA-DRFC, joan.decker@cea.fr)

0001 function [equil_out] = fitequil_yp(equil_in,mode,method,ngrid,nharm)
0002 % C3PO - Build vectorial tokamak magnetic equilibrium from its numerical counterpart
0003 %
0004 % Build vectorial tokamak magnetic equilibrium from its numerical counterpart
0005 %
0006 % Two operating modes :
0007 % 1) Spline-Fourier description for (psi,theta) magnetic equilibrium (HELENA for example)
0008 %
0009 %    1-D: use interp1.m MatLab built-in function (Hermite polynomials (pchip), spline, linear,  nearest,...)
0010 %    2-D: use Fourier series expansion for the poloidal dependence and the
0011 %    interp1.m for the Fourier series coefficients (Hermite polynomials (pchip), spline, linear,
0012 %    nearest,...)
0013 %
0014 % 2) Spline-Spline description for (x,y) magnetic equilibrium (EFIT for example, or for TORPEX case)
0015 %
0016 %    1-D: use interp1.m MatLab built-in function (Hermite polynomials (pchip), spline, linear, nearest,...)
0017 %    2-D: use interp1.m for the x-direction and the interp1.m again for the coefficients of the
0018 %    previous spline (Hermite polynomials (pchip), spline, bi-cubic, linear, nearest,...)
0019 %
0020 % By Yves Peysson (CEA-DRFC, yves.peysson@cea.fr) and Joan Decker (CEA-DRFC, joan.decker@cea.fr)
0021 %
0022 if nargin == 1, 
0023     mode = 1;%Spline-Fourier case
0024     method = 'spline';
0025     ngrid = 1001;
0026     nharm = [];
0027 end
0028 if nargin == 2, 
0029     method = 'spline';
0030     ngrid = 1001;
0031     nharm = [];
0032 end
0033 if nargin == 3, 
0034     ngrid = 1001;  
0035     nharm = [];
0036 end
0037 if nargin == 4, 
0038     nharm = [];%For Spline-Fourier case only
0039 end
0040 %
0041 equil_out.id = equil_in.id;
0042 equil_out.method = method;
0043 %
0044 equil_out.Rp = equil_in.Rp;
0045 equil_out.Zp = equil_in.Zp;
0046 %
0047 if mode == 1,%Spline-Fourier case: (psi,theta) magnetic equilibrium
0048     if isfield(equil_in,'pzni') && isfield(equil_in,'pzTi') && isfield(equil_in,'zZi') && isfield(equil_in,'pne') && isfield(equil_in,'pTe') 
0049         if ~isempty(equil_in.zZi),equil_out.zZi = equil_in.zZi;else equil_out.zZi = [0,0,0];end
0050         if ~isempty(equil_in.zmi),equil_out.zmi = equil_in.zmi;else equil_out.zmi = [1,1,1];end
0051     end
0052     %
0053     ap = equil_in.ptx(end,1);
0054     psi_apRp = equil_in.psi_apRp;
0055     psia_apRp = equil_in.psi_apRp(end);
0056     %
0057     sigmaIp =  sign(psia_apRp);%Sign of the plasma current (> 0 if Ip is clock-wise when the tokamak in viewed from the top)
0058     %
0059     %Resampling on a finer mesh close to psi_apRp = 0 before Hermite-Fourier calculations.
0060     %
0061     srho = linspace(0,1,ngrid).^2;
0062     %
0063     [rho] = psi2rho_jd(equil_in,psi_apRp/psia_apRp,method);
0064     %
0065     spsin = interp1(rho,psi_apRp/psia_apRp,srho,method);
0066     %
0067     if isfield(equil_in,'pzni') && isfield(equil_in,'pzTi') && isfield(equil_in,'zZi') && isfield(equil_in,'pne') && isfield(equil_in,'pTe') 
0068         if ~isempty(equil_in.pTe),sTe = interp1(rho,equil_in.pTe,srho,method);end
0069         if ~isempty(equil_in.pne),sne = interp1(rho,equil_in.pne,srho,method);end
0070         if ~isempty(equil_in.pzTi),szTi = interp1(rho,equil_in.pzTi',srho',method)';end
0071         if ~isempty(equil_in.pzni),szni = interp1(rho,equil_in.pzni',srho',method)';end
0072     end
0073     %
0074     sBx = interp1(rho,equil_in.ptBx,srho,method);%Bx field component
0075     sBy = interp1(rho,equil_in.ptBy,srho,method);%By field component
0076     sBz = interp1(rho,equil_in.ptBPHI,srho,method);%PHI field component
0077     %
0078     sBx(:,end) = sBx(:,1);
0079     sBy(:,end) = sBy(:,1);
0080     sBz(:,end) = sBz(:,1);
0081     %
0082     sx = interp1(rho,equil_in.ptx,srho,method);%x coordinate
0083     sy = interp1(rho,equil_in.pty,srho,method);%y coordinate
0084     %
0085     sx(:,end) = sx(:,1);
0086     sy(:,end) = sy(:,1);
0087     %
0088     sr = sqrt(sx.*sx + sy.*sy);%r
0089     %
0090     sBP = sqrt(sBx.*sBx + sBy.*sBy);%BP field
0091     sB = sqrt(sBx.*sBx +sBy.*sBy + sBz.*sBz);%B field
0092     %
0093     scalpha = (sx.*sBy - sy.*sBx)./sBP./sr;%cos(alpha)
0094     ssalpha = sigmaIp*(sx.*sBx + sy.*sBy)./sBP./sr;%sin(alpha)
0095     %
0096     scalpha(1,:) = scalpha(2,:);%to avoid singularity at r = 0
0097     ssalpha(1,:) = ssalpha(2,:);%to avoid singularity at r = 0
0098     %
0099     if isinf(equil_in.Rp),
0100         sUpsilon = 1;
0101     else
0102         sUpsilon = 1 + sx/equil_in.Rp;
0103     end
0104     %
0105     sdpsin = gradient(spsin(:));
0106     sdxstardpsin = gradient(sx(:,1))./sdpsin(:);
0107     %
0108     sgradrho = sUpsilon.*sBP.*(sdxstardpsin*ones(1,size(sx,2)))/abs(equil_in.psi_apRp(end));
0109     sgradrho(1,:) = sgradrho(2,:);
0110     sgradrho(:,end) = sgradrho(:,1);
0111     %
0112     % Calculation of the vectorial magnetic equilibrium
0113     %
0114     equil_out.psin_fit = interpequilpt_yp(method,spsin,srho);
0115     %
0116     sxstar_psin_fit = interpequilpt_yp(method,sx(:,1),spsin);%xstar = x(:,theta=0)
0117     sdxstardpsin_psin = ppval_yp(sxstar_psin_fit.pp_dfdrho,spsin);
0118     equil_out.dxstardpsin_fit = interpequilpt_yp(method,sdxstardpsin_psin,srho);
0119     %
0120     if isfield(equil_in,'pzni') && isfield(equil_in,'pzTi') && isfield(equil_in,'zZi') && isfield(equil_in,'pne') && isfield(equil_in,'pTe') 
0121         if ~isempty(equil_in.pTe),equil_out.Te_fit = interpequilpt_yp(method,sTe,srho);end
0122         if ~isempty(equil_in.pne),equil_out.ne_fit = interpequilpt_yp(method,sne,srho);end
0123         if ~isempty(equil_in.pzTi),equil_out.zTi_fit = interpequilpt_yp(method,szTi,srho);end
0124         if ~isempty(equil_in.pzni),equil_out.zni_fit = interpequilpt_yp(method,szni,srho);end
0125     end
0126     %
0127     if isempty(nharm) || isinf(nharm) || isnan(nharm) ,nharm = length(equil_in.theta) - 1;end
0128     %
0129     equil_out.Bx_fit = interpequilpt_yp(method,sBx,srho,equil_in.theta,nharm);%Bx field component
0130     equil_out.By_fit = interpequilpt_yp(method,sBy,srho,equil_in.theta,nharm);%By field component
0131     equil_out.Bz_fit = interpequilpt_yp(method,sBz,srho,equil_in.theta,nharm);%Bz field component
0132     equil_out.BP_fit = interpequilpt_yp(method,sBP,srho,equil_in.theta,nharm);%BP field
0133     equil_out.B_fit = interpequilpt_yp(method,sB,srho,equil_in.theta,nharm);%B field
0134     equil_out.x_fit = interpequilpt_yp(method,sx,srho,equil_in.theta,nharm);%x coordinate
0135     equil_out.y_fit = interpequilpt_yp(method,sy,srho,equil_in.theta,nharm);%y coordinate
0136     equil_out.r_fit = interpequilpt_yp(method,sr,srho,equil_in.theta,nharm);%y coordinate
0137     equil_out.calpha_fit = interpequilpt_yp(method,scalpha,srho,equil_in.theta,nharm);%cos(alpha)
0138     equil_out.salpha_fit = interpequilpt_yp(method,ssalpha,srho,equil_in.theta,nharm);%sin(alpha)
0139     equil_out.gradrho_fit = interpequilpt_yp(method,sgradrho,srho,equil_in.theta,nharm);%sin(alpha)
0140     %
0141     if psi_apRp(1) == 0,
0142         equil_out.ap = ap;
0143         equil_out.psia_apRp = psia_apRp;
0144     else
0145         error('equil_out.psi_apRp(1) must be equal to zero !');
0146     end
0147     %
0148 elseif mode == 2,%Spline-Spline description for (x,y) magnetic equilibrium
0149     %
0150     method = 'cubic';
0151     equil_out.method = method;
0152     %
0153     sxyx = linspace(min(equil_in.xyx),max(equil_in.xyx),ngrid);sxydx = sxyx(2) - sxyx(1);
0154     sxyy = linspace(max(equil_in.xyy),min(equil_in.xyy),ngrid);sxydy = sxyy(2) - sxyy(1);
0155     %
0156     [X,Y] = meshgrid(equil_in.xyx,equil_in.xyy);
0157     [sX,sY] = meshgrid(sxyx,sxyy);
0158     %
0159     sX_ij = sX(1:size(sX,1)-1,1:size(sX,2)-1);
0160     sX_ipj = sX(2:size(sX,1),1:size(sX,2)-1);
0161     sX_ijp = sX(1:size(sX,1)-1,2:size(sX,2));
0162     sX_ipjp = sX(1:size(sX,1)-1,2:size(sX,2));
0163     %
0164     sY_ij = sY(1:size(sY,1)-1,1:size(sY,2)-1);
0165     sY_ipj = sY(2:size(sY,1),1:size(sY,2)-1);
0166     sY_ijp = sY(1:size(sY,1)-1,2:size(sY,2));
0167     sY_ipjp = sY(1:size(sY,1)-1,2:size(sY,2));
0168     %
0169     equil_out.XX = [sX_ipj(:),sX_ipjp(:),sX_ijp(:),sX_ij(:)];
0170     equil_out.YY = [sY_ipj(:),sY_ipjp(:),sY_ijp(:),sY_ij(:)];
0171     %
0172     sxyne = interp2(X,Y,equil_in.xyne,sX,sY,method);
0173     [sxydnedx,sxydnedy] = gradient(sxyne,sxydx,sxydy);
0174     [sxyd2nedxdx,sxyd2nedxdy] = gradient(sxydnedx,sxydx,sxydy);
0175     [sxyd2nedydx,sxyd2nedydy] = gradient(sxydnedy,sxydx,sxydy);
0176     sxyd2nedxdy = (sxyd2nedxdy + sxyd2nedydx)/2;
0177     %
0178     sxyTe = interp2(X,Y,equil_in.xyTe,sX,sY,method);
0179     [sxydTedx,sxydTedy] = gradient(sxyTe,sxydx,sxydy);
0180     [sxyd2Tedxdx,sxyd2Tedxdy] = gradient(sxydTedx,sxydx,sxydy);
0181     [sxyd2Tedydx,sxyd2Tedydy] = gradient(sxydTedy,sxydx,sxydy);
0182     sxyd2Tedxdy = (sxyd2Tedxdy + sxyd2Tedydx)/2;
0183     %
0184     N = size(equil_in.xyne,1);
0185     sN = size(sxyne,1);
0186     %
0187     for iz = 1:length(equil_in.zZi),
0188         sxyzni(sN*(iz-1)+1:sN*iz,:) = interp2(X,Y,equil_in.xyzni(N*(iz-1)+1:N*iz,:),sX,sY,method);
0189         [sxyzdnidx(sN*(iz-1)+1:sN*iz,:),sxyzdnidy(sN*(iz-1)+1:sN*iz,:)] = gradient(sxyzni(sN*(iz-1)+1:sN*iz,:),sxydx,sxydy);
0190         [sxyzd2nidxdx(sN*(iz-1)+1:sN*iz,:),sxyzd2nidxdy(sN*(iz-1)+1:sN*iz,:)] = gradient(sxyzdnidx(sN*(iz-1)+1:sN*iz,:),sxydx,sxydy);
0191         [sxyzd2nidydx(sN*(iz-1)+1:sN*iz,:),sxyzd2nidydy(sN*(iz-1)+1:sN*iz,:)] = gradient(sxyzdnidy(sN*(iz-1)+1:sN*iz,:),sxydx,sxydy);
0192         sxyzd2nidxdy(sN*(iz-1)+1:sN*iz,:) = (sxyzd2nidxdy(sN*(iz-1)+1:sN*iz,:) + sxyzd2nidydx(sN*(iz-1)+1:sN*iz,:))/2;
0193         %
0194         sxyzTi(sN*(iz-1)+1:sN*iz,:) = interp2(X,Y,equil_in.xyzTi(N*(iz-1)+1:N*iz,:),sX,sY,method);
0195         [sxyzdTidx(sN*(iz-1)+1:sN*iz,:),sxyzdTidy(sN*(iz-1)+1:sN*iz,:)] = gradient(sxyzTi(sN*(iz-1)+1:sN*iz,:),sxydx,sxydy);
0196         [sxyzd2Tidxdx(sN*(iz-1)+1:sN*iz,:),sxyzd2Tidxdy(sN*(iz-1)+1:sN*iz,:)] = gradient(sxyzdTidx(sN*(iz-1)+1:sN*iz,:),sxydx,sxydy);
0197         [sxyzd2Tidydx(sN*(iz-1)+1:sN*iz,:),sxyzd2Tidydy(sN*(iz-1)+1:sN*iz,:)] = gradient(sxyzdTidy(sN*(iz-1)+1:sN*iz,:),sxydx,sxydy);
0198         sxyzd2Tidxdy(sN*(iz-1)+1:sN*iz,:) = (sxyzd2Tidxdy(sN*(iz-1)+1:sN*iz,:) + sxyzd2Tidydx(sN*(iz-1)+1:sN*iz,:))/2;
0199      end
0200     %
0201     sxyBx = interp2(X,Y,equil_in.xyBx,sX,sY,method);
0202     [sxydBxdx,sxydBxdy] = gradient(sxyBx,sxydx,sxydy);
0203     [sxyd2Bxdxdx,sxyd2Bxdxdy] = gradient(sxydBxdx,sxydx,sxydy);
0204     [sxyd2Bxdydx,sxyd2Bxdydy] = gradient(sxydBxdy,sxydx,sxydy);
0205     sxyd2Bxdxdy = (sxyd2Bxdxdy + sxyd2Bxdydx)/2;
0206     %
0207     sxyBy = interp2(X,Y,equil_in.xyBy,sX,sY,method);
0208     [sxydBydx,sxydBydy] = gradient(sxyBy,sxydx,sxydy);
0209     [sxyd2Bydxdx,sxyd2Bydxdy] = gradient(sxydBydx,sxydx,sxydy);
0210     [sxyd2Bydydx,sxyd2Bydydy] = gradient(sxydBydy,sxydx,sxydy);
0211     sxyd2Bydxdy = (sxyd2Bydxdy + sxyd2Bydydx)/2;
0212     %
0213     sxyBz = interp2(X,Y,equil_in.xyBz,sX,sY,method);
0214     [sxydBzdx,sxydBzdy] = gradient(sxyBz,sxydx,sxydy);
0215     [sxyd2Bzdxdx,sxyd2Bzdxdy] = gradient(sxydBzdx,sxydx,sxydy);
0216     [sxyd2Bzdydx,sxyd2Bzdydy] = gradient(sxydBzdy,sxydx,sxydy);
0217     sxyd2Bzdxdy = (sxyd2Bzdxdy + sxyd2Bzdydx)/2;
0218     %
0219     equil_out.ne_fit = interpequilxy_yp(sxyne,sxydnedx,sxydnedy,sxyd2nedxdy);
0220     equil_out.Te_fit = interpequilxy_yp(sxyTe,sxydTedx,sxydTedy,sxyd2Tedxdy);
0221     %
0222     equil_out.zni_fit = interpequilxy_yp(sxyzni,sxyzdnidx,sxyzdnidy,sxyzd2nidxdy);
0223     equil_out.zTi_fit = interpequilxy_yp(sxyzTi,sxyzdTidx,sxyzdTidy,sxyzd2Tidxdy);
0224     %
0225     equil_out.Bx_fit = interpequilxy_yp(sxyBx,sxydBxdx,sxydBxdy,sxyd2Bxdxdy);%Bx field component
0226     equil_out.By_fit = interpequilxy_yp(sxyBy,sxydBydx,sxydBydy,sxyd2Bydxdy);%By field component
0227     equil_out.Bz_fit = interpequilxy_yp(sxyBz,sxydBzdx,sxydBzdy,sxyd2Bzdxdy);%Bz (BPHI) field component
0228     %
0229 else
0230     error('Error: this option is not yet implemented in fitequil_yp.m');
0231 end
0232 %
0233 
0234 
0235