equilval_yp

 This function calculates all parameters of the axisymmetric magnetic equilibrium at any position in the plasma

   INPUT :

       - equil_fit: fitted axisymmetric magnetic equilibrium structure

   - if equil_fit corresponds to (rho,theta) coordinates

       - a: rho -> normalized minor radius (a.u.) [1,1]
       - b: theta -> poloidal angle (rad) [1,1]

   - if equil_fit corresponds to (x,y) coordinates

       - a: x -> radial coordinate (R - Rp) [1,1]
       - b: y -> vertical coordinate (Z - Zp) [1,1]

   OUTPUT :

   if input coordinates are (rho,theta,z):

       - equil_val.Te,equil_val.dTedrho,equil_val.d2Tedrho2: electron temperature and its derivatives [1,1]
       - equil_val.ne,equil_val.dnedrho,equil_val.d2nedrho2: electron desnity and its derivatives [1,1]
       - equil_val.zTi,equil_val.dzTidrho,equil_val.d2zTidrho2: ion temperatures and their derivatives [1,1]
       - equil_val.zni,equil_val.dznidrho,equil_val.d2znidrho2: ion densities and their derivatives [1,1]
       - equil_val.x,equil_val.dxdtheta,equil_val.d2xdtheta2,equil_val.dxdrho,equil_val.d2xdthetadrho,equil_val.d2xdrho2: horizontal coordinate and its derivatives [1,1]
       - equil_val.y,equil_val.dydtheta,equil_val.d2ydtheta2,equil_val.dydrho,equil_val.d2ydthetadrho,equil_val.d2ydrho2: vertical coordinate and its derivatives [1,1]
       - equil_val.r,equil_val.drdtheta,equil_val.d2rdtheta2,equil_val.drdrho,equil_val.d2rdthetadrho,equil_val.d2rdrho2: minor radius and its derivatives [1,1]
       - equil_val.salpha,equil_val.dsalphadtheta,equil_val.d2salphadtheta2,equil_val.dsalphadrho,equil_val.d2salphadthetadrho,equil_val.d2salphadrho2: sin(alpha) and its derivatives [1,1]
       - equil_val.calpha,equil_val.dcalphadtheta,equil_val.d2calphadtheta2,equil_val.dcalphadrho,equil_val.d2calphadthetadrho,equil_val.d2calphadrho2: cos(alpha) and its derivatives [1,1]
       - equil_val.Bx,equil_val.dBxdtheta,equil_val.d2Bxdtheta2,equil_val.dBxdrho,equil_val.d2Bxdthetadrho,equil_val.d2Bxdrho2: horizontal B-field and its derivatives
       - equil_val.By,equil_val.dBydtheta,equil_val.d2Bydtheta2,equil_val.dBydrho,equil_val.d2Bydthetadrho,equil_val.d2Bydrho2: vertical B-field and its derivatives [1,1]
       - equil_val.Bz,equil_val.dBzdtheta,equil_val.d2Bzdtheta2,equil_val.dBzdrho,equil_val.d2Bzdthetadrho,equil_val.d2Bzdrho2: toroidal or axial B-field and its derivatives (signed value) [1,1]
       - equil_val.BP,equil_val.dBPdtheta,equil_val.d2BPdtheta2,equil_val.dBPdrho,equil_val.d2BPdthetadrho,equil_val.d2BPdrho2: poloidal B-field and its derivatives [1,1]
       - equil_val.B,equil_val.dBdtheta,equil_val.d2Bdtheta2,equil_val.dBdrho,equil_val.d2Bdthetadrho,equil_val.d2Bdrho2: B-field and its derivatives [1,1]
       - equil_val.gradrho,equil_val.dgradrhodtheta,equil_val.d2gradrhodtheta2,equil_val.dgradrhodrho,equil_val.d2gradrhodthetadrho,equil_val.d2gradrhodrho2: normalized gradient of rho and its derivatives [1,1]
       - equil_val.Upsilon,equil_val.dUpsilondtheta,equil_val.d2Upsilondtheta2,equil_val.dUpsilondrho,equil_val.d2Upsilondthetadrho,equil_val.d2Upsilondrho2: toroidal form factor and its derivatives [1,1]

   if input coordinates are (x,y,z):


 by Y. Peysson DRFC-CEA <yves.peysson@cea.fr> and J. Decker DRFC-CEA <joan.decker@cea.fr>

0001 function [equil_val] = equilval_yp(equil_fit,a,b)
0002 %
0003 % This function calculates all parameters of the axisymmetric magnetic equilibrium at any position in the plasma
0004 %
0005 %   INPUT :
0006 %
0007 %       - equil_fit: fitted axisymmetric magnetic equilibrium structure
0008 %
0009 %   - if equil_fit corresponds to (rho,theta) coordinates
0010 %
0011 %       - a: rho -> normalized minor radius (a.u.) [1,1]
0012 %       - b: theta -> poloidal angle (rad) [1,1]
0013 %
0014 %   - if equil_fit corresponds to (x,y) coordinates
0015 %
0016 %       - a: x -> radial coordinate (R - Rp) [1,1]
0017 %       - b: y -> vertical coordinate (Z - Zp) [1,1]
0018 %
0019 %   OUTPUT :
0020 %
0021 %   if input coordinates are (rho,theta,z):
0022 %
0023 %       - equil_val.Te,equil_val.dTedrho,equil_val.d2Tedrho2: electron temperature and its derivatives [1,1]
0024 %       - equil_val.ne,equil_val.dnedrho,equil_val.d2nedrho2: electron desnity and its derivatives [1,1]
0025 %       - equil_val.zTi,equil_val.dzTidrho,equil_val.d2zTidrho2: ion temperatures and their derivatives [1,1]
0026 %       - equil_val.zni,equil_val.dznidrho,equil_val.d2znidrho2: ion densities and their derivatives [1,1]
0027 %       - equil_val.x,equil_val.dxdtheta,equil_val.d2xdtheta2,equil_val.dxdrho,equil_val.d2xdthetadrho,equil_val.d2xdrho2: horizontal coordinate and its derivatives [1,1]
0028 %       - equil_val.y,equil_val.dydtheta,equil_val.d2ydtheta2,equil_val.dydrho,equil_val.d2ydthetadrho,equil_val.d2ydrho2: vertical coordinate and its derivatives [1,1]
0029 %       - equil_val.r,equil_val.drdtheta,equil_val.d2rdtheta2,equil_val.drdrho,equil_val.d2rdthetadrho,equil_val.d2rdrho2: minor radius and its derivatives [1,1]
0030 %       - equil_val.salpha,equil_val.dsalphadtheta,equil_val.d2salphadtheta2,equil_val.dsalphadrho,equil_val.d2salphadthetadrho,equil_val.d2salphadrho2: sin(alpha) and its derivatives [1,1]
0031 %       - equil_val.calpha,equil_val.dcalphadtheta,equil_val.d2calphadtheta2,equil_val.dcalphadrho,equil_val.d2calphadthetadrho,equil_val.d2calphadrho2: cos(alpha) and its derivatives [1,1]
0032 %       - equil_val.Bx,equil_val.dBxdtheta,equil_val.d2Bxdtheta2,equil_val.dBxdrho,equil_val.d2Bxdthetadrho,equil_val.d2Bxdrho2: horizontal B-field and its derivatives
0033 %       - equil_val.By,equil_val.dBydtheta,equil_val.d2Bydtheta2,equil_val.dBydrho,equil_val.d2Bydthetadrho,equil_val.d2Bydrho2: vertical B-field and its derivatives [1,1]
0034 %       - equil_val.Bz,equil_val.dBzdtheta,equil_val.d2Bzdtheta2,equil_val.dBzdrho,equil_val.d2Bzdthetadrho,equil_val.d2Bzdrho2: toroidal or axial B-field and its derivatives (signed value) [1,1]
0035 %       - equil_val.BP,equil_val.dBPdtheta,equil_val.d2BPdtheta2,equil_val.dBPdrho,equil_val.d2BPdthetadrho,equil_val.d2BPdrho2: poloidal B-field and its derivatives [1,1]
0036 %       - equil_val.B,equil_val.dBdtheta,equil_val.d2Bdtheta2,equil_val.dBdrho,equil_val.d2Bdthetadrho,equil_val.d2Bdrho2: B-field and its derivatives [1,1]
0037 %       - equil_val.gradrho,equil_val.dgradrhodtheta,equil_val.d2gradrhodtheta2,equil_val.dgradrhodrho,equil_val.d2gradrhodthetadrho,equil_val.d2gradrhodrho2: normalized gradient of rho and its derivatives [1,1]
0038 %       - equil_val.Upsilon,equil_val.dUpsilondtheta,equil_val.d2Upsilondtheta2,equil_val.dUpsilondrho,equil_val.d2Upsilondthetadrho,equil_val.d2Upsilondrho2: toroidal form factor and its derivatives [1,1]
0039 %
0040 %   if input coordinates are (x,y,z):
0041 %
0042 %
0043 % by Y. Peysson DRFC-CEA <yves.peysson@cea.fr> and J. Decker DRFC-CEA <joan.decker@cea.fr>
0044 %
0045 if nargin < 4
0046     mode = 1;
0047 end
0048 %
0049 equil_val.Rp = equil_fit.Rp;
0050 equil_val.Zp = equil_fit.Zp;
0051 %
0052 if isfield(equil_fit.Bx_fit,'pp_a0'),% (rho,theta) coordinates
0053     %
0054     rho = a;
0055     theta = b;
0056     %
0057     equil_val.psin = ppval_yp(equil_fit.psin_fit.pp_f,rho);
0058     equil_val.dpsindrho = ppval_yp(equil_fit.psin_fit.pp_dfdrho,rho);
0059     equil_val.dxstardpsin = ppval_yp(equil_fit.dxstardpsin_fit.pp_f,rho);
0060     equil_val.ddxstardpsindrho = ppval_yp(equil_fit.dxstardpsin_fit.pp_dfdrho,rho);
0061     %
0062     if isfield(equil_fit,'Te_fit'),
0063       equil_val.Te = ppval_yp(equil_fit.Te_fit.pp_f,rho);
0064       equil_val.dTedrho = ppval_yp(equil_fit.Te_fit.pp_dfdrho,rho);
0065       equil_val.d2Tedrho2 = ppval_yp(equil_fit.Te_fit.pp_d2fdrho2,rho);
0066     else
0067       equil_val.Te = [];
0068       equil_val.dTedrho = [];
0069       equil_val.d2Tedrho2 = [];
0070     end
0071     %
0072     if isfield(equil_fit,'ne_fit'),
0073       equil_val.ne = ppval_yp(equil_fit.ne_fit.pp_f,rho);
0074       equil_val.dnedrho = ppval_yp(equil_fit.ne_fit.pp_dfdrho,rho);
0075       equil_val.d2nedrho2 = ppval_yp(equil_fit.ne_fit.pp_d2fdrho2,rho);
0076     else
0077       equil_val.ne = [];
0078       equil_val.dnedrho = [];
0079       equil_val.d2nedrho2 = [];  
0080     end
0081     %
0082     if isfield(equil_fit,'zTi_fit'),
0083       equil_val.zTi = ppval_yp(equil_fit.zTi_fit.pp_f,rho);
0084       equil_val.dzTidrho = ppval_yp(equil_fit.zTi_fit.pp_dfdrho,rho);
0085       equil_val.d2zTidrho2 = ppval_yp(equil_fit.zTi_fit.pp_d2fdrho2,rho);
0086     else
0087       equil_val.zTi = [];
0088       equil_val.dzTidrho = [];
0089       equil_val.d2zTidrho2 = [];
0090     end
0091     %
0092     if isfield(equil_fit,'zni_fit'),
0093       equil_val.zni = ppval_yp(equil_fit.zni_fit.pp_f,rho);
0094       equil_val.dznidrho = ppval_yp(equil_fit.zni_fit.pp_dfdrho,rho);
0095       equil_val.d2znidrho2 = ppval_yp(equil_fit.zni_fit.pp_d2fdrho2,rho);
0096     else
0097       equil_val.zni = [];
0098       equil_val.dznidrho = [];
0099       equil_val.d2znidrho2 = [];
0100     end
0101     %
0102     x_a0 = ppval_yp(equil_fit.x_fit.pp_a0,rho);
0103     x_an = ppval_yp(equil_fit.x_fit.pp_an,rho);
0104     x_bn = ppval_yp(equil_fit.x_fit.pp_bn,rho);
0105     x_da0drho = ppval_yp(equil_fit.x_fit.pp_da0drho,rho);
0106     x_dandrho = ppval_yp(equil_fit.x_fit.pp_dandrho,rho);
0107     x_dbndrho = ppval_yp(equil_fit.x_fit.pp_dbndrho,rho);
0108     if isfield(equil_fit.x_fit,'pp_d2a0drho2'),
0109         x_d2a0drho2 = ppval_yp(equil_fit.x_fit.pp_d2a0drho2,rho);
0110         x_d2andrho2 = ppval_yp(equil_fit.x_fit.pp_d2andrho2,rho);
0111         x_d2bndrho2 = ppval_yp(equil_fit.x_fit.pp_d2bndrho2,rho);
0112     else
0113         x_d2a0drho2 = zeros(size(x_da0drho));
0114         x_d2andrho2 = zeros(size(x_dandrho));
0115         x_d2bndrho2 = zeros(size(x_dbndrho));
0116     end
0117     %
0118     [equil_val.x,equil_val.dxdtheta,equil_val.d2xdtheta2,equil_val.dxdrho,equil_val.d2xdthetadrho,equil_val.d2xdrho2] = calcval_yp(equil_fit,theta,x_a0,x_an,x_bn,x_da0drho,x_dandrho,x_dbndrho,x_d2a0drho2,x_d2andrho2,x_d2bndrho2);
0119     %
0120     y_a0 = ppval_yp(equil_fit.y_fit.pp_a0,rho);
0121     y_an = ppval_yp(equil_fit.y_fit.pp_an,rho);
0122     y_bn = ppval_yp(equil_fit.y_fit.pp_bn,rho);
0123     y_da0drho = ppval_yp(equil_fit.y_fit.pp_da0drho,rho);
0124     y_dandrho = ppval_yp(equil_fit.y_fit.pp_dandrho,rho);
0125     y_dbndrho = ppval_yp(equil_fit.y_fit.pp_dbndrho,rho);
0126     if isfield(equil_fit.y_fit,'pp_d2a0drho2');
0127         y_d2a0drho2 = ppval_yp(equil_fit.y_fit.pp_d2a0drho2,rho);
0128         y_d2andrho2 = ppval_yp(equil_fit.y_fit.pp_d2andrho2,rho);
0129         y_d2bndrho2 = ppval_yp(equil_fit.y_fit.pp_d2bndrho2,rho);
0130     else
0131         y_d2a0drho2 = zeros(size(y_da0drho));
0132         y_d2andrho2 = zeros(size(y_dandrho));
0133         y_d2bndrho2 =zeros(size(y_dbndrho));
0134     end
0135     %
0136     [equil_val.y,equil_val.dydtheta,equil_val.d2ydtheta2,equil_val.dydrho,equil_val.d2ydthetadrho,equil_val.d2ydrho2] = calcval_yp(equil_fit,theta,y_a0,y_an,y_bn,y_da0drho,y_dandrho,y_dbndrho,y_d2a0drho2,y_d2andrho2,y_d2bndrho2);
0137     %
0138     r_a0 = ppval_yp(equil_fit.r_fit.pp_a0,rho);
0139     r_an = ppval_yp(equil_fit.r_fit.pp_an,rho);
0140     r_bn = ppval_yp(equil_fit.r_fit.pp_bn,rho);
0141     r_da0drho = ppval_yp(equil_fit.r_fit.pp_da0drho,rho);
0142     r_dandrho = ppval_yp(equil_fit.r_fit.pp_dandrho,rho);
0143     r_dbndrho = ppval_yp(equil_fit.r_fit.pp_dbndrho,rho);
0144     if isfield(equil_fit.r_fit,'pp_d2a0drho2');
0145         r_d2a0drho2 = ppval_yp(equil_fit.r_fit.pp_d2a0drho2,rho);
0146         r_d2andrho2 = ppval_yp(equil_fit.r_fit.pp_d2andrho2,rho);
0147         r_d2bndrho2 = ppval_yp(equil_fit.r_fit.pp_d2bndrho2,rho);
0148     else
0149         r_d2a0drho2 = zeros(size(r_da0drho));
0150         r_d2andrho2 = zeros(size(r_dandrho));
0151         r_d2bndrho2 = zeros(size(r_dbndrho)); 
0152     end
0153     %
0154     [equil_val.r,equil_val.drdtheta,equil_val.d2rdtheta2,equil_val.drdrho,equil_val.d2rdthetadrho,equil_val.d2rdrho2] = calcval_yp(equil_fit,theta,r_a0,r_an,r_bn,r_da0drho,r_dandrho,r_dbndrho,r_d2a0drho2,r_d2andrho2,r_d2bndrho2);
0155     %
0156     Bx_a0 = ppval_yp(equil_fit.Bx_fit.pp_a0,rho);
0157     Bx_an = ppval_yp(equil_fit.Bx_fit.pp_an,rho);
0158     Bx_bn = ppval_yp(equil_fit.Bx_fit.pp_bn,rho);
0159     Bx_da0drho = ppval_yp(equil_fit.Bx_fit.pp_da0drho,rho);
0160     Bx_dandrho = ppval_yp(equil_fit.Bx_fit.pp_dandrho,rho);
0161     Bx_dbndrho = ppval_yp(equil_fit.Bx_fit.pp_dbndrho,rho);
0162     if isfield(equil_fit.Bx_fit,'pp_d2a0drho2');
0163         Bx_d2a0drho2 = ppval_yp(equil_fit.Bx_fit.pp_d2a0drho2,rho);
0164         Bx_d2andrho2 = ppval_yp(equil_fit.Bx_fit.pp_d2andrho2,rho);
0165         Bx_d2bndrho2 = ppval_yp(equil_fit.Bx_fit.pp_d2bndrho2,rho);
0166     else
0167         Bx_d2a0drho2 = zeros(size(Bx_da0drho));
0168         Bx_d2andrho2 = zeros(size(Bx_dandrho));
0169         Bx_d2bndrho2 = zeros(size(Bx_dbndrho)); 
0170     end
0171     %
0172     [equil_val.Bx,equil_val.dBxdtheta,equil_val.d2Bxdtheta2,equil_val.dBxdrho,equil_val.d2Bxdthetadrho,equil_val.d2Bxdrho2] = calcval_yp(equil_fit,theta,Bx_a0,Bx_an,Bx_bn,Bx_da0drho,Bx_dandrho,Bx_dbndrho,Bx_d2a0drho2,Bx_d2andrho2,Bx_d2bndrho2);
0173     %
0174     By_a0 = ppval_yp(equil_fit.By_fit.pp_a0,rho);
0175     By_an = ppval_yp(equil_fit.By_fit.pp_an,rho);
0176     By_bn = ppval_yp(equil_fit.By_fit.pp_bn,rho);
0177     By_da0drho = ppval_yp(equil_fit.By_fit.pp_da0drho,rho);
0178     By_dandrho = ppval_yp(equil_fit.By_fit.pp_dandrho,rho);
0179     By_dbndrho = ppval_yp(equil_fit.By_fit.pp_dbndrho,rho);
0180     if isfield(equil_fit.By_fit,'pp_d2a0drho2');
0181         By_d2a0drho2 = ppval_yp(equil_fit.By_fit.pp_d2a0drho2,rho);
0182         By_d2andrho2 = ppval_yp(equil_fit.By_fit.pp_d2andrho2,rho);
0183         By_d2bndrho2 = ppval_yp(equil_fit.By_fit.pp_d2bndrho2,rho);
0184     else
0185         By_d2a0drho2 = zeros(size(By_da0drho));
0186         By_d2andrho2 = zeros(size(By_dandrho));
0187         By_d2bndrho2 = zeros(size(By_dbndrho)); 
0188     end
0189     %
0190     [equil_val.By,equil_val.dBydtheta,equil_val.d2Bydtheta2,equil_val.dBydrho,equil_val.d2Bydthetadrho,equil_val.d2Bydrho2] = calcval_yp(equil_fit,theta,By_a0,By_an,By_bn,By_da0drho,By_dandrho,By_dbndrho,By_d2a0drho2,By_d2andrho2,By_d2bndrho2);
0191     %
0192     Bz_a0 = ppval_yp(equil_fit.Bz_fit.pp_a0,rho);
0193     Bz_an = ppval_yp(equil_fit.Bz_fit.pp_an,rho);
0194     Bz_bn = ppval_yp(equil_fit.Bz_fit.pp_bn,rho);
0195     Bz_da0drho = ppval_yp(equil_fit.Bz_fit.pp_da0drho,rho);
0196     Bz_dandrho = ppval_yp(equil_fit.Bz_fit.pp_dandrho,rho);
0197     Bz_dbndrho = ppval_yp(equil_fit.Bz_fit.pp_dbndrho,rho);
0198     if isfield(equil_fit.Bz_fit,'pp_d2a0drho2');
0199         Bz_d2a0drho2 = ppval_yp(equil_fit.Bz_fit.pp_d2a0drho2,rho);
0200         Bz_d2andrho2 = ppval_yp(equil_fit.Bz_fit.pp_d2andrho2,rho);
0201         Bz_d2bndrho2 = ppval_yp(equil_fit.Bz_fit.pp_d2bndrho2,rho);
0202     else
0203         Bz_d2a0drho2 = zeros(size(Bz_da0drho));
0204         Bz_d2andrho2 = zeros(size(Bz_dandrho));
0205         Bz_d2bndrho2 = zeros(size(Bz_dbndrho)); 
0206     end
0207     %
0208     [equil_val.Bz,equil_val.dBzdtheta,equil_val.d2Bzdtheta2,equil_val.dBzdrho,equil_val.d2Bzdthetadrho,equil_val.d2Bzdrho2] = calcval_yp(equil_fit,theta,Bz_a0,Bz_an,Bz_bn,Bz_da0drho,Bz_dandrho,Bz_dbndrho,Bz_d2a0drho2,Bz_d2andrho2,Bz_d2bndrho2);
0209     %
0210     BP_a0 = ppval_yp(equil_fit.BP_fit.pp_a0,rho);
0211     BP_an = ppval_yp(equil_fit.BP_fit.pp_an,rho);
0212     BP_bn = ppval_yp(equil_fit.BP_fit.pp_bn,rho);
0213     BP_da0drho = ppval_yp(equil_fit.BP_fit.pp_da0drho,rho);
0214     BP_dandrho = ppval_yp(equil_fit.BP_fit.pp_dandrho,rho);
0215     BP_dbndrho = ppval_yp(equil_fit.BP_fit.pp_dbndrho,rho);
0216     if isfield(equil_fit.BP_fit,'pp_d2a0drho2');
0217         BP_d2a0drho2 = ppval_yp(equil_fit.BP_fit.pp_d2a0drho2,rho);
0218         BP_d2andrho2 = ppval_yp(equil_fit.BP_fit.pp_d2andrho2,rho);
0219         BP_d2bndrho2 = ppval_yp(equil_fit.BP_fit.pp_d2bndrho2,rho);
0220     else
0221         BP_d2a0drho2 = zeros(size(BP_da0drho));
0222         BP_d2andrho2 = zeros(size(BP_dandrho));
0223         BP_d2bndrho2 = zeros(size(BP_dbndrho)); 
0224     end
0225     %
0226     [equil_val.BP,equil_val.dBPdtheta,equil_val.d2BPdtheta2,equil_val.dBPdrho,equil_val.d2BPdthetadrho,equil_val.d2BPdrho2] = calcval_yp(equil_fit,theta,BP_a0,BP_an,BP_bn,BP_da0drho,BP_dandrho,BP_dbndrho,BP_d2a0drho2,BP_d2andrho2,BP_d2bndrho2);
0227     %
0228     B_a0 = ppval_yp(equil_fit.B_fit.pp_a0,rho);
0229     B_an = ppval_yp(equil_fit.B_fit.pp_an,rho);
0230     B_bn = ppval_yp(equil_fit.B_fit.pp_bn,rho);
0231     B_da0drho = ppval_yp(equil_fit.B_fit.pp_da0drho,rho);
0232     B_dandrho = ppval_yp(equil_fit.B_fit.pp_dandrho,rho);
0233     B_dbndrho = ppval_yp(equil_fit.B_fit.pp_dbndrho,rho);
0234     if isfield(equil_fit.B_fit,'pp_d2a0drho2');
0235         B_d2a0drho2 = ppval_yp(equil_fit.B_fit.pp_d2a0drho2,rho);
0236         B_d2andrho2 = ppval_yp(equil_fit.B_fit.pp_d2andrho2,rho);
0237         B_d2bndrho2 = ppval_yp(equil_fit.B_fit.pp_d2bndrho2,rho);
0238     else
0239         B_d2a0drho2 = zeros(size(B_da0drho));
0240         B_d2andrho2 = zeros(size(B_dandrho));
0241         B_d2bndrho2 = zeros(size(B_dbndrho)); 
0242     end
0243     %
0244     [equil_val.B,equil_val.dBdtheta,equil_val.d2Bdtheta2,equil_val.dBdrho,equil_val.d2Bdthetadrho,equil_val.d2Bdrho2] = calcval_yp(equil_fit,theta,B_a0,B_an,B_bn,B_da0drho,B_dandrho,B_dbndrho,B_d2a0drho2,B_d2andrho2,B_d2bndrho2);
0245     %
0246     calpha_a0 = ppval_yp(equil_fit.calpha_fit.pp_a0,rho);
0247     calpha_an = ppval_yp(equil_fit.calpha_fit.pp_an,rho);
0248     calpha_bn = ppval_yp(equil_fit.calpha_fit.pp_bn,rho);
0249     calpha_da0drho = ppval_yp(equil_fit.calpha_fit.pp_da0drho,rho);
0250     calpha_dandrho = ppval_yp(equil_fit.calpha_fit.pp_dandrho,rho);
0251     calpha_dbndrho = ppval_yp(equil_fit.calpha_fit.pp_dbndrho,rho);
0252     if isfield(equil_fit.calpha_fit,'pp_d2a0drho2');
0253         calpha_d2a0drho2 = ppval_yp(equil_fit.calpha_fit.pp_d2a0drho2,rho);
0254         calpha_d2andrho2 = ppval_yp(equil_fit.calpha_fit.pp_d2andrho2,rho);
0255         calpha_d2bndrho2 = ppval_yp(equil_fit.calpha_fit.pp_d2bndrho2,rho);
0256     else
0257         calpha_d2a0drho2 = zeros(size(calpha_da0drho));
0258         calpha_d2andrho2 = zeros(size(calpha_dandrho));
0259         calpha_d2bndrho2 = zeros(size(calpha_dbndrho)); 
0260     end
0261     %
0262     [equil_val.calpha,equil_val.dcalphadtheta,equil_val.d2calphadtheta2,equil_val.dcalphadrho,equil_val.d2calphadthetadrho,equil_val.d2calphadrho2] = calcval_yp(equil_fit,theta,calpha_a0,calpha_an,calpha_bn,calpha_da0drho,calpha_dandrho,calpha_dbndrho,calpha_d2a0drho2,calpha_d2andrho2,calpha_d2bndrho2);
0263     %
0264     salpha_a0 = ppval_yp(equil_fit.salpha_fit.pp_a0,rho);
0265     salpha_an = ppval_yp(equil_fit.salpha_fit.pp_an,rho);
0266     salpha_bn = ppval_yp(equil_fit.salpha_fit.pp_bn,rho);
0267     salpha_da0drho = ppval_yp(equil_fit.salpha_fit.pp_da0drho,rho);
0268     salpha_dandrho = ppval_yp(equil_fit.salpha_fit.pp_dandrho,rho);
0269     salpha_dbndrho = ppval_yp(equil_fit.salpha_fit.pp_dbndrho,rho);
0270     if isfield(equil_fit.salpha_fit,'pp_d2a0drho2');
0271         salpha_d2a0drho2 = ppval_yp(equil_fit.salpha_fit.pp_d2a0drho2,rho);
0272         salpha_d2andrho2 = ppval_yp(equil_fit.salpha_fit.pp_d2andrho2,rho);
0273         salpha_d2bndrho2 = ppval_yp(equil_fit.salpha_fit.pp_d2bndrho2,rho);
0274     else
0275         salpha_d2a0drho2 = zeros(size(salpha_da0drho));
0276         salpha_d2andrho2 = zeros(size(salpha_dandrho));
0277         salpha_d2bndrho2 = zeros(size(salpha_dbndrho)); 
0278     end
0279     %
0280     [equil_val.salpha,equil_val.dsalphadtheta,equil_val.d2salphadtheta2,equil_val.dsalphadrho,equil_val.d2salphadthetadrho,equil_val.d2salphadrho2] = calcval_yp(equil_fit,theta,salpha_a0,salpha_an,salpha_bn,salpha_da0drho,salpha_dandrho,salpha_dbndrho,salpha_d2a0drho2,salpha_d2andrho2,salpha_d2bndrho2);
0281     %
0282     gradrho_a0 = ppval_yp(equil_fit.gradrho_fit.pp_a0,rho);
0283     gradrho_an = ppval_yp(equil_fit.gradrho_fit.pp_an,rho);
0284     gradrho_bn = ppval_yp(equil_fit.gradrho_fit.pp_bn,rho);
0285     gradrho_da0drho = ppval_yp(equil_fit.gradrho_fit.pp_da0drho,rho);
0286     gradrho_dandrho = ppval_yp(equil_fit.gradrho_fit.pp_dandrho,rho);
0287     gradrho_dbndrho = ppval_yp(equil_fit.gradrho_fit.pp_dbndrho,rho);
0288     if isfield(equil_fit.gradrho_fit,'pp_d2a0drho2');
0289         gradrho_d2a0drho2 = ppval_yp(equil_fit.gradrho_fit.pp_d2a0drho2,rho);
0290         gradrho_d2andrho2 = ppval_yp(equil_fit.gradrho_fit.pp_d2andrho2,rho);
0291         gradrho_d2bndrho2 = ppval_yp(equil_fit.gradrho_fit.pp_d2bndrho2,rho);
0292     else
0293         gradrho_d2a0drho2 = zeros(size(gradrho_da0drho));
0294         gradrho_d2andrho2 = zeros(size(gradrho_dandrho));
0295         gradrho_d2bndrho2 = zeros(size(gradrho_dbndrho)); 
0296     end
0297     %
0298     [equil_val.gradrho,equil_val.dgradrhodtheta,equil_val.d2gradrhodtheta2,equil_val.dgradrhodrho,equil_val.d2gradrhodthetadrho,equil_val.d2gradrhodrho2] = calcval_yp(equil_fit,theta,gradrho_a0,gradrho_an,gradrho_bn,gradrho_da0drho,gradrho_dandrho,gradrho_dbndrho,gradrho_d2a0drho2,gradrho_d2andrho2,gradrho_d2bndrho2);
0299     %
0300     if isinf(equil_val.Rp),
0301         equil_val.Upsilon = 1;
0302         equil_val.dUpsilondrho = 0;
0303         equil_val.dUpsilondtheta = 0;
0304         equil_val.d2Upsilondtheta2 = 0;
0305         equil_val.d2Upsilondthetadrho = 0;
0306         equil_val.d2Upsilondrho2 = 0;
0307     else
0308 
0309         equil_val.Upsilon = 1 + equil_val.x/equil_val.Rp;
0310         equil_val.dUpsilondrho = equil_val.dxdrho/equil_val.Rp;
0311         equil_val.dUpsilondtheta = equil_val.dxdtheta/equil_val.Rp;
0312         equil_val.d2Upsilondtheta2 = equil_val.d2xdtheta2/equil_val.Rp;
0313         equil_val.d2Upsilondthetadrho = equil_val.d2xdthetadrho/equil_val.Rp;
0314         equil_val.d2Upsilondrho2 = equil_val.d2xdrho2/equil_val.Rp;
0315     end
0316     %
0317 elseif isfield(equil_fit.Bx_fit,'CC'),% (x,y,phi) coordinates
0318     %
0319     if size(a) == size(b),
0320         x = a(:).';
0321         y = b(:).';
0322     else
0323         error('ERROR: input coordinates must have the same number of dimensions');
0324     end
0325     %
0326 
0327     keyboard
0328     [k,l] = find((equil_fit.XX(:,1)*ones(1,length(x)) < ones(length(equil_fit.XX(:,1)),1)*x) & (ones(length(equil_fit.XX(:,1)),1)*x < equil_fit.XX(:,2)*ones(1,length(x))) & (equil_fit.YY(:,1)*ones(1,length(y)) < ones(length(equil_fit.YY(:,1)),1)*y) & (ones(length(equil_fit.YY(:,4)),1)*y < equil_fit.YY(:,4)*ones(1,length(y))));
0329     %
0330     t = ((x.' - equil_fit.XX(k,1))./(equil_fit.XX(k,2) - equil_fit.XX(k,1))).';
0331     u = ((y.' - equil_fit.YY(k,1))./(equil_fit.YY(k,4) - equil_fit.YY(k,1))).';
0332     %
0333     l = 1;
0334     for ii = 0:3,
0335         for jj = 0:3,
0336             ut(l,:) = t.^ii.*u.^jj;
0337             dutdt(l,:) = ii*t.^(ii-1).*u.^jj;
0338             dutdu(l,:) = jj*t.^ii.*u.^(jj-1);
0339             l = l+1;
0340         end
0341     end
0342     %
0343     equil_val.ne = equil_fit.ne_fit.CC(k,:)*ut(:);
0344     equil_val.dnedx = equil_fit.ne_fit.CC(k,:)*dutdt(:);
0345     equil_val.dnedy = equil_fit.ne_fit.CC(k,:)*dutdu(:);
0346     %
0347 
0348 else
0349     error('ERROR: the format of the structure equil_fit is not recognized')
0350 end