interpequilxy_yp

 C3PO: Bicubic interpolation of the magnetic equilibrium

 Build interpolated toroidal MHD magnetic equilibrium from its numerical counterpart
 Bicubic interpolation (x = R-Rp, y = Z-Zp). 

 INPUT

    - f_in: 2-D (x,y) function to interpolate [nz*nx,nx]
    - dfdx_in: 2-D (x,y) 1st order derivative of the function to interpolate along x direction [nz*nx,nx]
    - dfdy_in: 2-D (x,y) 1st order derivative of the function to interpolate along y direction [nz*nx,nx]
    - d2fdxdy_in: 2-D (x,y) 2nd order cross-derivative of the function to interpolate along x and y directions [nz*nx,nx]

 OUTPUT

    - Zbci: structure of bicubic interpolation

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

0001 function [Zbci] = interpequilxy_yp(f_in,dfdx_in,dfdy_in,d2fdxdy_in)
0002 % C3PO: Bicubic interpolation of the magnetic equilibrium
0003 %
0004 % Build interpolated toroidal MHD magnetic equilibrium from its numerical counterpart
0005 % Bicubic interpolation (x = R-Rp, y = Z-Zp).
0006 %
0007 % INPUT
0008 %
0009 %    - f_in: 2-D (x,y) function to interpolate [nz*nx,nx]
0010 %    - dfdx_in: 2-D (x,y) 1st order derivative of the function to interpolate along x direction [nz*nx,nx]
0011 %    - dfdy_in: 2-D (x,y) 1st order derivative of the function to interpolate along y direction [nz*nx,nx]
0012 %    - d2fdxdy_in: 2-D (x,y) 2nd order cross-derivative of the function to interpolate along x and y directions [nz*nx,nx]
0013 %
0014 % OUTPUT
0015 %
0016 %    - Zbci: structure of bicubic interpolation
0017 %
0018 % By Yves Peysson (CEA-DRFC, yves.peysson@cea.fr) and Joan Decker (CEA-DRFC, joan.decker@epfl.ch)
0019 %
0020 if nargin < 4,
0021     error('ERROR: not enough input arguments in interpequilxy_yp');
0022 end
0023 %
0024 wt_d = [1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0;
0025 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0;
0026 -3, 0, 0, 3, 0, 0, 0, 0,-2, 0, 0,-1, 0, 0, 0, 0;
0027 2, 0, 0,-2, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0;
0028 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0;
0029 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0;
0030 0, 0, 0, 0,-3, 0, 0, 3, 0, 0, 0, 0,-2, 0, 0,-1;
0031 0, 0, 0, 0, 2, 0, 0,-2, 0, 0, 0, 0, 1, 0, 0, 1;
0032 -3, 3, 0, 0,-2,-1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0;
0033 0, 0, 0, 0, 0, 0, 0, 0,-3, 3, 0, 0,-2,-1, 0, 0;
0034 9,-9, 9,-9, 6, 3,-3,-6, 6,-6,-3, 3, 4, 2, 1, 2;
0035 -6, 6,-6, 6,-4,-2, 2, 4,-3, 3, 3,-3,-2,-1,-1,-2;
0036 2,-2, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0;
0037 0, 0, 0, 0, 0, 0, 0, 0, 2,-2, 0, 0, 1, 1, 0, 0;
0038 -6, 6,-6, 6,-3,-3, 3, 3,-4, 4, 2,-2,-2,-2,-1,-1;
0039 4,-4, 4,-4, 2, 2,-2,-2, 2,-2,-2, 2, 1, 1, 1, 1];%weights (see note on redmine)
0040 %
0041 nwt_d = size(wt_d,1);
0042 %
0043 for jj = 1:size(f_in,1)/size(f_in,2),
0044 %
0045     f = f_in((jj-1)*size(f_in,2)+1:jj*size(f_in,2),:);
0046     dfdx = dfdx_in((jj-1)*size(f_in,2)+1:jj*size(f_in,2),:);
0047     dfdy = dfdy_in((jj-1)*size(f_in,2)+1:jj*size(f_in,2),:);
0048     d2fdxdy = d2fdxdy_in((jj-1)*size(f_in,2)+1:jj*size(f_in,2),:);
0049     %
0050     F_ij = f(1:size(f,1)-1,1:size(f,2)-1);
0051     F_ipj = f(2:size(f,1),1:size(f,2)-1);
0052     F_ijp = f(1:size(f,1)-1,2:size(f,2));
0053     F_ipjp = f(1:size(f,1)-1,2:size(f,2));
0054     %
0055     dFdx_ij = dfdx(1:size(f,1)-1,1:size(f,2)-1);
0056     dFdx_ipj = dfdx(2:size(f,1),1:size(f,2)-1);
0057     dFdx_ijp = dfdx(1:size(f,1)-1,2:size(f,2));
0058     dFdx_ipjp = dfdx(1:size(f,1)-1,2:size(f,2));
0059     %
0060     dFdy_ij = dfdy(1:size(f,1)-1,1:size(f,2)-1);
0061     dFdy_ipj = dfdy(2:size(f,1),1:size(f,2)-1);
0062     dFdy_ijp = dfdy(1:size(f,1)-1,2:size(f,2));
0063     dFdy_ipjp = dfdy(1:size(f,1)-1,2:size(f,2));
0064     %
0065     d2Fdxdy_ij = d2fdxdy(1:size(f,1)-1,1:size(f,2)-1);
0066     d2Fdxdy_ipj = d2fdxdy(2:size(f,1),1:size(f,2)-1);
0067     d2Fdxdy_ijp = d2fdxdy(1:size(f,1)-1,2:size(f,2));
0068     d2Fdxdy_ipjp = d2fdxdy(1:size(f,1)-1,2:size(f,2));
0069     %
0070     NN = size(F_ij,1)*size(F_ij,2);
0071     %
0072     FF = [F_ipj(:),F_ipjp(:),F_ijp(:),F_ij(:),dFdx_ipj(:),dFdx_ipjp(:),dFdx_ijp(:),dFdx_ij(:),dFdy_ipj(:),dFdy_ipjp(:),dFdy_ijp(:),dFdy_ij(:),d2Fdxdy_ipj(:),d2Fdxdy_ipjp(:),d2Fdxdy_ijp(:),d2Fdxdy_ij(:)].';
0073     FF = FF(:);
0074     %
0075     clear F_ij F_ipj F_ijp F_ipjp dFdx_ij dFdx_ipj dFdx_ijp dFdx_ipjp dFdy_ij dFdy_ipj dFdy_ijp dFdy_ipjp d2Fdxdy_ij d2Fdxdy_ipj d2Fdxdy_ijp d2Fdxdy_ipjp
0076     %
0077     %M_test = kron(speye(NN),sparse(wt_d));%not efficient when NN is very large
0078     %CC_test = reshape(M_test*FF,nwt_d,NN).';;%Create bicubic interpolation coefficients
0079     %
0080     CC = zeros(NN*nwt_d,1);
0081     NM0 = min(NN,1000);%1000 gives minium time to calculate CC usually.
0082     M0 = sparse(NM0*nwt_d,NM0*nwt_d);
0083     M0 = kron(speye(NM0),sparse(wt_d));
0084     %
0085     for ii = 1:floor(NN*nwt_d/size(M0,1))
0086         CC(size(M0,1)*(ii-1)+1:size(M0,1)*ii,1) = M0*FF(size(M0,1)*(ii-1)+1:size(M0,1)*ii,1);%Create bicubic interpolation coefficients
0087     end
0088     %
0089     NM1 = NN - floor(NN*nwt_d/size(M0,1))*NM0;
0090     M1 = sparse(NM1*nwt_d,NM1*nwt_d);
0091     M1 = kron(speye(NM1),sparse(wt_d));
0092     %
0093     CC(NN*nwt_d-size(M1,1)+1:NN*nwt_d,1) = M1*FF(NN*nwt_d-size(M1,1)+1:NN*nwt_d,1);%Create last bicubic interpolation coefficients
0094     %
0095     Zbci.CC(:,:,jj) = reshape(CC,nwt_d,NN).';
0096     %
0097 end