interpequilpt_yp

PURPOSE ^

SYNOPSIS ^

function [Zpp_f] = interpequilpt_yp(method,f,rho,theta,nharm)

DESCRIPTION ^

 Build interpolated toroidal MHD magnetic equilibrium from its numerical counterpart
 Piecewise polynomials (1D, radial) or piecewise polynomials + Fourier series (2D, radial + poloidal)

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

 INPUT
    - method: interpolation method -> Hermite polynomial(pchip), spline, nearest, linear, bicubic, ...) [nchar,1]
    - f: 1-D (rho) or 2-D (rho,theta) function to interpolate [nrho,1] or [nrho,mtheta]
    - rho: normalized minor radius [nrho,1]
    - theta: poloidal flux grid coordinate [1,ntheta]
    - nharm: number of harmonics for the Fourier series expansion (<ntheta-1) [1,1]

 OUTPUT

    - Zpp_f: structure of piecewise polynomials (1-D) and Fourier series coefficients (2-D)

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

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function [Zpp_f] = interpequilpt_yp(method,f,rho,theta,nharm)
0002 %
0003 % Build interpolated toroidal MHD magnetic equilibrium from its numerical counterpart
0004 % Piecewise polynomials (1D, radial) or piecewise polynomials + Fourier series (2D, radial + poloidal)
0005 %
0006 % 1-D: use interp1 MatLab built-in function (Hermite polynomial (pchip) spline, bi-cubic, linear, nearest,...)
0007 % 2-D: use Fourier series expansion for the poloidal dependence and the interp1 for the Fourier series coefficients (spline, bi-cubic, linear, nearest,...)
0008 %
0009 % INPUT
0010 %    - method: interpolation method -> Hermite polynomial(pchip), spline, nearest, linear, bicubic, ...) [nchar,1]
0011 %    - f: 1-D (rho) or 2-D (rho,theta) function to interpolate [nrho,1] or [nrho,mtheta]
0012 %    - rho: normalized minor radius [nrho,1]
0013 %    - theta: poloidal flux grid coordinate [1,ntheta]
0014 %    - nharm: number of harmonics for the Fourier series expansion (<ntheta-1) [1,1]
0015 %
0016 % OUTPUT
0017 %
0018 %    - Zpp_f: structure of piecewise polynomials (1-D) and Fourier series coefficients (2-D)
0019 %
0020 % By Yves Peysson (CEA-DRFC, yves.peysson@cea.fr) and Joan Decker (CEA-DRFC, joan.decker@cea.fr)
0021 %
0022 if nargin == 3,
0023     %
0024     if isvector(f),f = f(:)';end
0025     Zpp_f.pp_f = ppinterp_yp(rho(:),f',size(f'),method);%Create piecewise polynomials
0026     %
0027     if ~isfield(Zpp_f.pp_f,'orient'),
0028         Zpp_f.pp_f.orient = 'first';
0029     end
0030     %
0031     Zpp_f.pp_dfdrho = Zpp_f.pp_f;%1st order derivative along rho
0032     Zpp_f.pp_d2fdrho2 = Zpp_f.pp_f;%2nd order derivative along rho
0033     %
0034     if Zpp_f.pp_f.order == 2,%Linear
0035         Zpp_f.pp_dfdrho.coefs = [Zpp_f.pp_f.coefs(:,1)];%1st order partial derivative along rho
0036         Zpp_f.pp_d2fdrho2.coefs = zeros(size([Zpp_f.pp_f.coefs(:,1)]));%2nd order partial derivative along rho
0037     elseif Zpp_f.pp_f.order == 3,%Quadratic
0038         Zpp_f.pp_dfdrho.coefs = [2*Zpp_f.pp_f.coefs(:,1),Zpp_f.pp_f.coefs(:,2)];%1st order partial derivative along rho
0039         Zpp_f.pp_d2fdrho2.coefs = [2*Zpp_f.pp_f.coefs(:,1)];%2nd order partial derivative along rho
0040     elseif Zpp_f.pp_f.order == 4,%cubic spline or Hermite polynomial (pchip option)
0041         Zpp_f.pp_dfdrho.coefs = [3*Zpp_f.pp_f.coefs(:,1),2*Zpp_f.pp_f.coefs(:,2),Zpp_f.pp_f.coefs(:,3)];%1st order partial derivative along rho
0042         Zpp_f.pp_d2fdrho2.coefs = [6*Zpp_f.pp_f.coefs(:,1),2*Zpp_f.pp_f.coefs(:,2)];%2nd order partial derivative along rho
0043     end 
0044     %
0045     Zpp_f.pp_dfdrho.order = Zpp_f.pp_f.order - 1;%1st order derivative along rho
0046     Zpp_f.pp_d2fdrho2.order = Zpp_f.pp_f.order - 2;%2nd order derivative along rho
0047     %
0048 elseif nargin >= 4,
0049     if min(size(f)) == 1, error('Matrix dimensions do not agree !');end
0050     rho = rho(:);
0051     theta = theta(:);
0052     %
0053     if nargin == 4,nharm = length(theta)-1;end
0054     %
0055     if nargin == 5 & nharm > length(theta)-1,
0056         %nharm = length(theta)-1;
0057         %info_dke_yp(2,['WARNING: the number of harmonics may not exceed : ',int2str(length(theta)-1)]);
0058     end
0059     %
0060     if nargin == 5 & isempty(nharm),
0061         nharm = length(theta)-1;
0062     end
0063     %
0064     FScoef = zeros(size(f,1),2*nharm + 1);
0065     for ii = 1:size(f,1),
0066         FSfun = fscreate(theta,f(ii,:),nharm,'foh');%Fourier series coefficients calculations
0067         FScoef(ii,:) = FSfun('coef');   
0068     end
0069     %
0070     Zpp_f.n = [1:nharm]';%Number of harmonics
0071     %
0072     a0 = FScoef(:,nharm+1);%f(rho,theta) = a0(rho) + sum_n(an(rho)*cos(n*theta) + bn(rho)*sin(n*theta);
0073     an = 2*real(FScoef(:,[nharm+2:end]));
0074     bn = -2*imag(FScoef(:,[nharm+2:end]));
0075     %
0076     Zpp_f.pp_a0 = ppinterp_yp(rho,a0,size(a0),method);%Interpolation of Fourier series coefficients on a fine rho grid
0077     Zpp_f.pp_an = ppinterp_yp(rho,an,size(an),method);%Interpolation of Fourier series coefficients on a fine rho grid
0078     Zpp_f.pp_bn = ppinterp_yp(rho,bn,size(bn),method);%Interpolation of Fourier series coefficients on a fine rho grid
0079     %
0080     if ~isfield(Zpp_f.pp_a0,'orient'),
0081         Zpp_f.pp_a0.orient = 'first';
0082         Zpp_f.pp_an.orient = 'first';
0083         Zpp_f.pp_bn.orient = 'first';
0084     end
0085     %
0086     Zpp_f.pp_da0drho = Zpp_f.pp_a0;
0087     Zpp_f.pp_dandrho = Zpp_f.pp_an;
0088     Zpp_f.pp_dbndrho = Zpp_f.pp_bn;
0089     Zpp_f.pp_d2a0drho2 = Zpp_f.pp_a0;
0090     Zpp_f.pp_d2andrho2 = Zpp_f.pp_an;
0091     Zpp_f.pp_d2bndrho2 = Zpp_f.pp_bn;
0092     %
0093     if Zpp_f.pp_a0.order == 2,%Linear
0094         Zpp_f.pp_da0drho.coefs = [Zpp_f.pp_a0.coefs(:,1)];%1st order partial derivative along rho
0095         Zpp_f.pp_dandrho.coefs = [Zpp_f.pp_an.coefs(:,1)];%1st order partial derivative along rho
0096         Zpp_f.pp_dbndrho.coefs = [Zpp_f.pp_bn.coefs(:,1)];%1st order partial derivative along rho
0097         Zpp_f.pp_d2a0drho2.coefs = zeros(size([Zpp_f.pp_a0.coefs(:,1)]));%2nd order partial derivative along rho
0098         Zpp_f.pp_d2andrho2.coefs = zeros(size([Zpp_f.pp_an.coefs(:,1)]));%2nd order partial derivative along rho
0099         Zpp_f.pp_d2bndrho2.coefs = zeros(size([Zpp_f.pp_bn.coefs(:,1)]));%2nd order partial derivative along rho
0100     elseif Zpp_f.pp_a0.order == 3,%Quadratic
0101         Zpp_f.pp_da0drho.coefs = [2*Zpp_f.pp_a0.coefs(:,1),Zpp_f.pp_a0.coefs(:,2)];%1st order partial derivative along rho
0102         Zpp_f.pp_dandrho.coefs = [2*Zpp_f.pp_an.coefs(:,1),Zpp_f.pp_an.coefs(:,2)];%1st order partial derivative along rho
0103         Zpp_f.pp_dbndrho.coefs = [2*Zpp_f.pp_bn.coefs(:,1),Zpp_f.pp_bn.coefs(:,2)];%1st order partial derivative along rho
0104         Zpp_f.pp_d2a0drho2.coefs = [2*Zpp_f.pp_a0.coefs(:,1)];%2nd order partial derivative along rho
0105         Zpp_f.pp_d2andrho2.coefs = [2*Zpp_f.pp_an.coefs(:,1)];%2nd order partial derivative along rho
0106         Zpp_f.pp_d2bndrho2.coefs = [2*Zpp_f.pp_bn.coefs(:,1)];%2nd order partial derivative along rho
0107     elseif Zpp_f.pp_a0.order == 4,%cubic spline or Hermite polynomials (option pchip)
0108         Zpp_f.pp_da0drho.coefs = [3*Zpp_f.pp_a0.coefs(:,1),2*Zpp_f.pp_a0.coefs(:,2),Zpp_f.pp_a0.coefs(:,3)];%1st order partial derivative along rho
0109         Zpp_f.pp_dandrho.coefs = [3*Zpp_f.pp_an.coefs(:,1),2*Zpp_f.pp_an.coefs(:,2),Zpp_f.pp_an.coefs(:,3)];%1st order partial derivative along rho
0110         Zpp_f.pp_dbndrho.coefs = [3*Zpp_f.pp_bn.coefs(:,1),2*Zpp_f.pp_bn.coefs(:,2),Zpp_f.pp_bn.coefs(:,3)];%1st order partial derivative along rho
0111         Zpp_f.pp_d2a0drho2.coefs = [6*Zpp_f.pp_a0.coefs(:,1),2*Zpp_f.pp_a0.coefs(:,2)];%2nd order partial derivative along rho
0112         Zpp_f.pp_d2andrho2.coefs = [6*Zpp_f.pp_an.coefs(:,1),2*Zpp_f.pp_an.coefs(:,2)];%2nd order partial derivative along rho
0113         Zpp_f.pp_d2bndrho2.coefs = [6*Zpp_f.pp_bn.coefs(:,1),2*Zpp_f.pp_bn.coefs(:,2)];%2nd order partial derivative along rho
0114     end
0115     %
0116     Zpp_f.pp_da0drho.order = Zpp_f.pp_a0.order - 1;%1st order partial derivative along rho
0117     Zpp_f.pp_dandrho.order = Zpp_f.pp_an.order - 1;%1st order partial derivative along rho
0118     Zpp_f.pp_dbndrho.order = Zpp_f.pp_bn.order - 1;%1st order partial derivative along rho
0119     Zpp_f.pp_d2a0drho2.order = Zpp_f.pp_a0.order - 2;%2nd order partial derivative along rho
0120     Zpp_f.pp_d2andrho2.order = Zpp_f.pp_an.order - 2;%2nd order partial derivative along rho
0121     Zpp_f.pp_d2bndrho2.order = Zpp_f.pp_bn.order - 2;%2nd order partial derivative along rho
0122 else
0123     error('Coordinates are missing !');
0124 end
Community support and wiki are available on Redmine. Last update: 18-Apr-2019.