testfitequil_yp

 Test the interpolated magnetic equilibrium

 INPUT
    - equil_in: numerical magnetic equilibrium structure
    - equil_fit: interpolated magnetic equilibrium structure
    - rho: normalized minor radius [1,1]
      (default = 0.25)
    - theta: poloidal angle (rad) [1,1]
      (default = 0.2)
    - ns: number of species [1,1]
      (default = 2)

 OUTPUT: (size(rho)=size(theta)=(1,1) only)
    - R: Major axis position [1,1]
    - Z: Vertical axis position [1,1]

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

0001 function [x,y] = testfitequil_yp(equil_in,equil_fit,rho,theta,ns)
0002 %
0003 % Test the interpolated magnetic equilibrium
0004 %
0005 % INPUT
0006 %    - equil_in: numerical magnetic equilibrium structure
0007 %    - equil_fit: interpolated magnetic equilibrium structure
0008 %    - rho: normalized minor radius [1,1]
0009 %      (default = 0.25)
0010 %    - theta: poloidal angle (rad) [1,1]
0011 %      (default = 0.2)
0012 %    - ns: number of species [1,1]
0013 %      (default = 2)
0014 %
0015 % OUTPUT: (size(rho)=size(theta)=(1,1) only)
0016 %    - R: Major axis position [1,1]
0017 %    - Z: Vertical axis position [1,1]
0018 %
0019 % By Yves Peysson (CEA-DRFC, yves.peysson@cea.fr) and Joan Decker (CEA-DRFC, joan.decker@cea.fr)
0020 %
0021 close all
0022 %
0023 flagcalc = 0;
0024 %
0025 if nargin < 3,% Output for comparison with C raytracing_yp.c
0026    rho = 0.25;
0027    theta = 0.2;
0028    ns = 2;
0029    flagcalc = 1;
0030 end
0031 if nargin < 4,% Output for comparison with C raytracing_yp.c
0032    theta = 0.2;
0033    ns = 2;
0034 end
0035 if nargin < 5,
0036    ns = 2;
0037 end
0038 %
0039 x = [];
0040 y = [];
0041 %
0042 prec = 1e-3;
0043 %
0044 sigmaIp = sign(equil_in.psi_apRp(end));
0045 %
0046 psin = equil_in.psi_apRp/equil_in.psi_apRp(end);
0047 [prho] = psi2rho_jd(equil_in,psin,'spline');
0048 dpsin = gradient(psin(:));
0049 d2psin = gradient(dpsin(:));
0050 dprho = gradient(prho(:));
0051 dtheta = gradient(equil_in.theta);
0052 %
0053 pTe = equil_in.pTe(:);
0054 pne = equil_in.pne(:);
0055 pzTi = equil_in.pzTi';
0056 pzni = equil_in.pzni';
0057 %
0058 Bx = equil_in.ptBx + eps;
0059 By = equil_in.ptBy + eps;
0060 Bz = equil_in.ptBPHI;
0061 %
0062 Rp = equil_in.Rp;
0063 Zp = equil_in.Zp;
0064 x = equil_in.ptx;
0065 y = equil_in.pty;
0066 %
0067 BP = sqrt(Bx.^2 + By.^2);
0068 B = sqrt(Bx.^2 + By.^2 + Bz.^2);
0069 r = sqrt(x.^2 + y.^2)+eps;
0070 calpha = (x.*By - y.*Bx)./BP./r;%cos(alpha)
0071 calpha(1,:) = calpha(2,:);%cos(alpha)
0072 salpha = sigmaIp*(x.*Bx + y.*By)./BP./r;%cos(alpha)
0073 salpha(1,:) = salpha(2,:);%cos(alpha)
0074 %
0075 % Numerical calculation of derivatives
0076 %
0077 [drdtheta,drdrho] = gradient(r);%r
0078 drdtheta = drdtheta./(ones(length(dprho),1)*dtheta);%1st order derivative
0079 drdrho = drdrho./(dprho(:)*ones(1,length(dtheta)));%1st order derivative
0080 [d2rdtheta2,d2rdthetadrho] = gradient(drdtheta);
0081 [d2rdrhodtheta,d2rdrho2] = gradient(drdrho);
0082 d2rdtheta2 = d2rdtheta2./(ones(length(dprho),1)*dtheta);%2nd order derivative
0083 d2rdthetadrho = d2rdthetadrho./(dprho(:)*ones(1,length(dtheta)));%2nd order derivative
0084 d2rdrhodtheta = d2rdrhodtheta./(ones(length(dprho),1)*dtheta);%2nd order derivative
0085 d2rdrho2 = d2rdrho2./(dprho(:)*ones(1,length(dtheta)));%2nd order derivative
0086 %
0087 [dxdtheta,dxdrho] = gradient(x);%Xx
0088 dxdtheta = dxdtheta./(ones(length(dprho),1)*dtheta);
0089 dxdrho = dxdrho./(dprho(:)*ones(1,length(dtheta)));
0090 [d2xdtheta2,d2xdthetadrho] = gradient(dxdtheta);
0091 [d2xdrhodtheta,d2xdrho2] = gradient(dxdrho);
0092 d2xdtheta2 = d2xdtheta2./(ones(length(dprho),1)*dtheta);%2nd order derivative
0093 d2xdthetadrho = d2xdthetadrho./(dprho(:)*ones(1,length(dtheta)));%2nd order derivative
0094 d2xdrhodtheta = d2xdrhodtheta./(ones(length(dprho),1)*dtheta);%2nd order derivative
0095 d2xdrho2 = d2xdrho2./(dprho(:)*ones(1,length(dtheta)));%2nd order derivative
0096 %
0097 [dydtheta,dydrho] = gradient(y);%Xy
0098 dydtheta = dydtheta./(ones(length(dprho),1)*dtheta);
0099 dydrho = dydrho./(dprho(:)*ones(1,length(dtheta)));
0100 [d2ydtheta2,d2ydthetadrho] = gradient(dydtheta);
0101 [d2ydrhodtheta,d2ydrho2] = gradient(dydrho);
0102 d2ydtheta2 = d2ydtheta2./(ones(length(dprho),1)*dtheta);%2nd order derivative
0103 d2ydthetadrho = d2ydthetadrho./(dprho(:)*ones(1,length(dtheta)));%2nd order derivative
0104 d2xdrhodtheta = d2ydrhodtheta./(ones(length(dprho),1)*dtheta);%2nd order derivative
0105 d2ydrho2 = d2ydrho2./(dprho(:)*ones(1,length(dtheta)));%2nd order derivative
0106 %
0107 [dcalphadtheta,dcalphadrho] = gradient(calpha);%
0108 dcalphadtheta = dcalphadtheta./(ones(length(dprho),1)*dtheta);
0109 dcalphadrho = dcalphadrho./(dprho(:)*ones(1,length(dtheta)));
0110 [d2calphadtheta2,d2calphadthetadrho] = gradient(dcalphadtheta);
0111 [d2calphadrhodtheta,d2calphadrho2] = gradient(dcalphadrho);
0112 d2calphadtheta2 = d2calphadtheta2./(ones(length(dprho),1)*dtheta);%2nd order derivative
0113 d2calphadthetadrho = d2calphadthetadrho./(dprho(:)*ones(1,length(dtheta)));%2nd order derivative
0114 d2calphadrhodtheta = d2calphadrhodtheta./(ones(length(dprho),1)*dtheta);%2nd order derivative
0115 d2calphadrho2 = d2calphadrho2./(dprho(:)*ones(1,length(dtheta)));%2nd order derivative
0116 %
0117 [dsalphadtheta,dsalphadrho] = gradient(salpha);%
0118 dsalphadtheta = dsalphadtheta./(ones(length(dprho),1)*dtheta);
0119 dsalphadrho = dsalphadrho./(dprho(:)*ones(1,length(dtheta)));
0120 [d2salphadtheta2,d2salphadthetadrho] = gradient(dsalphadtheta);
0121 [d2salphadrhodtheta,d2salphadrho2] = gradient(dsalphadrho);
0122 d2salphadtheta2 = d2salphadtheta2./(ones(length(dprho),1)*dtheta);%2nd order derivative
0123 d2salphadthetadrho = d2salphadthetadrho./(dprho(:)*ones(1,length(dtheta)));%2nd order derivative
0124 d2salphadrhodtheta = d2salphadrhodtheta./(ones(length(dprho),1)*dtheta);%2nd order derivative
0125 d2salphadrho2 = d2salphadrho2./(dprho(:)*ones(1,length(dtheta)));%2nd order derivative
0126 %
0127 [dBzdtheta,dBzdrho] = gradient(Bz);%Toroidal magnetic field
0128 dBzdtheta = dBzdtheta./(ones(length(dprho),1)*dtheta);
0129 dBzdrho = dBzdrho./(dprho(:)*ones(1,length(dtheta)));
0130 [d2Bzdtheta2,d2Bzdthetadrho] = gradient(dBzdtheta);
0131 [d2Bzdrhodtheta,d2Bzdrho2] = gradient(dBzdrho);
0132 d2Bzdtheta2 = d2Bzdtheta2./(ones(length(dprho),1)*dtheta);%2nd order derivative
0133 d2Bzdthetadrho = d2Bzdthetadrho./(dprho(:)*ones(1,length(dtheta)));%2nd order derivative
0134 d2Bzdrhodtheta = d2Bzdrhodtheta./(ones(length(dprho),1)*dtheta);%2nd order derivative
0135 d2Bzdrho2 = d2Bzdrho2./(dprho(:)*ones(1,length(dtheta)));%2nd order derivative
0136 %
0137 [dBPdtheta,dBPdrho] = gradient(BP);%Poloidal magnetic field
0138 dBPdtheta = dBPdtheta./(ones(length(dprho),1)*dtheta);
0139 dBPdrho = dBPdrho./(dprho(:)*ones(1,length(dtheta)));
0140 [d2BPdtheta2,d2BPdthetadrho] = gradient(dBPdtheta);
0141 [d2BPdrhodtheta,d2BPdrho2] = gradient(dBPdrho);
0142 d2BPdtheta2 = d2BPdtheta2./(ones(length(dprho),1)*dtheta);%2nd order derivative
0143 d2BPdthetadrho = d2BPdthetadrho./(dprho(:)*ones(1,length(dtheta)));%2nd order derivative
0144 d2BPdrhodtheta = d2BPdrhodtheta./(ones(length(dprho),1)*dtheta);%2nd order derivative
0145 d2BPdrho2 = d2BPdrho2./(dprho(:)*ones(1,length(dtheta)));%2nd order derivative
0146 %
0147 [dBdtheta,dBdrho] = gradient(B);%Total magnetic field
0148 dBdtheta = dBdtheta./(ones(length(dprho),1)*dtheta);
0149 dBdrho = dBdrho./(dprho(:)*ones(1,length(dtheta)));
0150 [d2Bdtheta2,d2Bdthetadrho] = gradient(dBdtheta);
0151 [d2Bdrhodtheta,d2Bdrho2] = gradient(dBdrho);
0152 d2Bdtheta2 = d2Bdtheta2./(ones(length(dprho),1)*dtheta);%2nd order derivative
0153 d2Bdthetadrho = d2Bdthetadrho./(dprho(:)*ones(1,length(dtheta)));%2nd order derivative
0154 d2Bdrhodtheta = d2Bdrhodtheta./(ones(length(dprho),1)*dtheta);%2nd order derivative
0155 d2Bdrho2 = d2Bdrho2./(dprho(:)*ones(1,length(dtheta)));%2nd order derivative
0156 %
0157 dpsindrho = dpsin./dprho;
0158 d2psindrho2 = d2psin./dprho./dprho;
0159 %
0160 if ~isempty(pTe),
0161     pdTe = gradient(pTe);
0162     pd2Te = gradient(pdTe);
0163     dTedrho = pdTe./dprho(:);
0164     d2Tedrho2 = pd2Te./dprho(:)./dprho(:);
0165 end
0166 if ~isempty(pne),
0167     pdne = gradient(pne);
0168     pd2ne = gradient(pdne);
0169     dnedrho = pdne./dprho(:);
0170     d2nedrho2 = pd2ne./dprho(:)./dprho(:);
0171 end
0172 if ~isempty(equil_in.zZi),
0173    for i = 1:length(equil_in.zZi)
0174       pdzTi = gradient(pzTi(:,i));
0175       pd2zTi = gradient(pdzTi);       
0176       pdzni = gradient(pzni(:,i));
0177       pd2zni = gradient(pdzni);       
0178       %
0179       zdTidrho(:,i) = pdzTi./dprho(:);
0180       zdnidrho(:,i) = pdzni./dprho(:);
0181       %
0182       zd2Tidrho2(:,i) = pd2zTi./dprho(:)./dprho(:);
0183       zd2nidrho2(:,i) = pd2zni./dprho(:)./dprho(:);
0184    end
0185 end
0186 %
0187 dxstardpsin = gradient(x(:,1))./dpsin(:);
0188 %
0189 if isinf(Rp),
0190     Upsilon = 1;
0191 else
0192     Upsilon = 1 + x/Rp;
0193 end
0194 %
0195 gradrho = Upsilon.*BP.*(dxstardpsin*ones(1,size(x,2)))/abs(equil_fit.psia_apRp);
0196 %
0197 [dgradrhodtheta,dgradrhodrho] = gradient(gradrho);%rho gradient
0198 dgradrhodtheta = dgradrhodtheta./(ones(length(dprho),1)*dtheta);
0199 dgradrhodrho = dgradrhodrho./(dprho(:)*ones(1,length(dtheta)));
0200 [d2gradrhodtheta2,d2gradrhodthetadrho] = gradient(dgradrhodtheta);
0201 [d2gradrhodrhodtheta,d2gradrhodrho2] = gradient(dgradrhodrho);
0202 d2gradrhodtheta2 = d2gradrhodtheta2./(ones(length(dprho),1)*dtheta);%2nd order derivative
0203 d2gradrhodthetadrho = d2gradrhodthetadrho./(dprho(:)*ones(1,length(dtheta)));%2nd order derivative
0204 d2gradrhodrhodtheta = d2gradrhodrhodtheta./(ones(length(dprho),1)*dtheta);%2nd order derivative
0205 d2gradrhodrho2 = d2gradrhodrho2./(dprho(:)*ones(1,length(dtheta)));%2nd order derivative
0206 %
0207 % Calculation of derivatives with Matlab build-in interpolation using polynomials (spline...)
0208 %
0209 flagcalc = 0
0210 if flagcalc == 0
0211     %
0212    close all;
0213 %
0214    nrho_fit = 500;
0215    ntheta_fit= length(equil_in.theta);
0216 %
0217    rho_fit = linspace(prho(1),prho(end),nrho_fit)';
0218    rho_fit(1) = prec;
0219    theta_fit = equil_in.theta';
0220 %
0221    psin_fit = ppval(equil_fit.psin_fit.pp_f,rho_fit);psin_fit = psin_fit(:);
0222    psin_fit1 = ppval_yp(equil_fit.psin_fit.pp_f,rho_fit);psin_fit1 = psin_fit1(:);
0223    dpsindrho_fit = ppval(equil_fit.psin_fit.pp_dfdrho,rho_fit);
0224    dpsindrho_fit1 = ppval_yp(equil_fit.psin_fit.pp_dfdrho,rho_fit);
0225    d2psindrho2_fit = ppval(equil_fit.psin_fit.pp_d2fdrho2,rho_fit);
0226    d2psindrho2_fit1 = ppval_yp(equil_fit.psin_fit.pp_d2fdrho2,rho_fit);
0227 %
0228    if isfield(equil_fit,'Te_fit'),
0229       Te_fit = ppval(equil_fit.Te_fit.pp_f,rho_fit);Te_fit = Te_fit(:);
0230       Te_fit1 = ppval_yp(equil_fit.Te_fit.pp_f,rho_fit);Te_fit1 = Te_fit1(:);
0231       dTedrho_fit = ppval(equil_fit.Te_fit.pp_dfdrho,rho_fit);
0232       dTedrho_fit1 = ppval_yp(equil_fit.Te_fit.pp_dfdrho,rho_fit);
0233       d2Tedrho2_fit = ppval(equil_fit.Te_fit.pp_d2fdrho2,rho_fit);
0234       d2Tedrho2_fit1 = ppval_yp(equil_fit.Te_fit.pp_d2fdrho2,rho_fit);
0235    end
0236 %
0237    if isfield(equil_fit,'ne_fit'),
0238       ne_fit = ppval(equil_fit.ne_fit.pp_f,rho_fit);ne_fit = ne_fit(:);
0239       ne_fit1 = ppval_yp(equil_fit.ne_fit.pp_f,rho_fit);ne_fit1 = ne_fit1(:);
0240       dnedrho_fit = ppval(equil_fit.ne_fit.pp_dfdrho,rho_fit);
0241       dnedrho_fit1 = ppval_yp(equil_fit.ne_fit.pp_dfdrho,rho_fit);
0242       d2nedrho2_fit = ppval(equil_fit.ne_fit.pp_d2fdrho2,rho_fit);
0243       d2nedrho2_fit1 = ppval_yp(equil_fit.ne_fit.pp_d2fdrho2,rho_fit);
0244    end
0245 %
0246    if isfield(equil_fit,'zTi_fit'),
0247       zTi_fit = ppval(equil_fit.zTi_fit.pp_f,rho_fit);
0248       zTi_fit1 = ppval_yp(equil_fit.zTi_fit.pp_f,rho_fit);
0249       zdTidrho_fit = ppval(equil_fit.zTi_fit.pp_dfdrho,rho_fit);
0250       zdTidrho_fit1 = ppval_yp(equil_fit.zTi_fit.pp_dfdrho,rho_fit);
0251       zd2Tidrho2_fit = ppval(equil_fit.zTi_fit.pp_d2fdrho2,rho_fit);
0252       zd2Tidrho2_fit1 = ppval_yp(equil_fit.zTi_fit.pp_d2fdrho2,rho_fit);
0253    end
0254 %
0255    if isfield(equil_fit,'zni_fit'),
0256       zni_fit = ppval(equil_fit.zni_fit.pp_f,rho_fit);
0257       zni_fit1 = ppval_yp(equil_fit.zni_fit.pp_f,rho_fit);
0258       zdnidrho_fit = ppval(equil_fit.zni_fit.pp_dfdrho,rho_fit);
0259       zdnidrho_fit1 = ppval_yp(equil_fit.zni_fit.pp_dfdrho,rho_fit);
0260       zd2nidrho2_fit = ppval(equil_fit.zni_fit.pp_d2fdrho2,rho_fit);
0261       zd2nidrho2_fit1 = ppval_yp(equil_fit.zni_fit.pp_d2fdrho2,rho_fit);
0262    end
0263 %
0264 % Calculation of the coordinates and magnetic equilibrium quantity
0265 %
0266    nr = 5; 
0267    rho_fit1 = prho(1:nr:end);%New normalized radial grid for display 2-D magnetic field
0268    rho_fit1(1) = prec;
0269 %
0270    dxstardpsin_fit = ppval(equil_fit.dxstardpsin_fit.pp_f,rho_fit1);dxstardpsin_fit = dxstardpsin_fit(:);
0271    dxstardpsin_fit1 = ppval_yp(equil_fit.dxstardpsin_fit.pp_f,rho_fit1);dxstardpsin_fit1 = dxstardpsin_fit1(:); 
0272    ddxstardpsindrho_fit = ppval(equil_fit.dxstardpsin_fit.pp_dfdrho,rho_fit1);ddxstardpsindrho_fit = ddxstardpsindrho_fit(:);
0273    ddxstardpsindrho_fit1 = ppval_yp(equil_fit.dxstardpsin_fit.pp_dfdrho,rho_fit1);ddxstardpsindrho_fit1 = ddxstardpsindrho_fit1(:);
0274    d2dxstardpsindrho2_fit = ppval(equil_fit.dxstardpsin_fit.pp_d2fdrho2,rho_fit1);d2dxstardpsindrho2_fit = d2dxstardpsindrho2_fit(:);
0275    d2dxstardpsindrho2_fit1 = ppval_yp(equil_fit.dxstardpsin_fit.pp_d2fdrho2,rho_fit1);d2dxstardpsindrho2_fit1 = d2dxstardpsindrho2_fit1(:);
0276 %
0277    x_a0_fit = ppval(equil_fit.x_fit.pp_a0,rho_fit1);
0278    x_an_fit = ppval(equil_fit.x_fit.pp_an,rho_fit1);
0279    x_bn_fit = ppval(equil_fit.x_fit.pp_bn,rho_fit1);
0280    x_da0drho_fit = ppval(equil_fit.x_fit.pp_da0drho,rho_fit1);
0281    x_dandrho_fit = ppval(equil_fit.x_fit.pp_dandrho,rho_fit1);
0282    x_dbndrho_fit = ppval(equil_fit.x_fit.pp_dbndrho,rho_fit1);
0283    x_d2a0drho2_fit = ppval(equil_fit.x_fit.pp_d2a0drho2,rho_fit1);
0284    x_d2andrho2_fit = ppval(equil_fit.x_fit.pp_d2andrho2,rho_fit1);
0285    x_d2bndrho2_fit = ppval(equil_fit.x_fit.pp_d2bndrho2,rho_fit1);
0286 %
0287    y_a0_fit = ppval(equil_fit.y_fit.pp_a0,rho_fit1);
0288    y_an_fit = ppval(equil_fit.y_fit.pp_an,rho_fit1);
0289    y_bn_fit = ppval(equil_fit.y_fit.pp_bn,rho_fit1);
0290    y_da0drho_fit = ppval(equil_fit.y_fit.pp_da0drho,rho_fit1);
0291    y_dandrho_fit = ppval(equil_fit.y_fit.pp_dandrho,rho_fit1);
0292    y_dbndrho_fit = ppval(equil_fit.y_fit.pp_dbndrho,rho_fit1);
0293    y_d2a0drho2_fit = ppval(equil_fit.y_fit.pp_d2a0drho2,rho_fit1);
0294    y_d2andrho2_fit = ppval(equil_fit.y_fit.pp_d2andrho2,rho_fit1);
0295    y_d2bndrho2_fit = ppval(equil_fit.y_fit.pp_d2bndrho2,rho_fit1);
0296 %
0297    Bx_a0_fit = ppval(equil_fit.Bx_fit.pp_a0,rho_fit1);
0298    Bx_an_fit = ppval(equil_fit.Bx_fit.pp_an,rho_fit1);
0299    Bx_bn_fit = ppval(equil_fit.Bx_fit.pp_bn,rho_fit1);
0300    Bx_da0drho_fit = ppval(equil_fit.Bx_fit.pp_da0drho,rho_fit1);
0301    Bx_dandrho_fit = ppval(equil_fit.Bx_fit.pp_dandrho,rho_fit1);
0302    Bx_dbndrho_fit = ppval(equil_fit.Bx_fit.pp_dbndrho,rho_fit1);
0303    Bx_d2a0drho2_fit = ppval(equil_fit.Bx_fit.pp_d2a0drho2,rho_fit1);
0304    Bx_d2andrho2_fit = ppval(equil_fit.Bx_fit.pp_d2andrho2,rho_fit1);
0305    Bx_d2bndrho2_fit = ppval(equil_fit.Bx_fit.pp_d2bndrho2,rho_fit1);
0306 %
0307    By_a0_fit = ppval(equil_fit.By_fit.pp_a0,rho_fit1);
0308    By_an_fit = ppval(equil_fit.By_fit.pp_an,rho_fit1);
0309    By_bn_fit = ppval(equil_fit.By_fit.pp_bn,rho_fit1);
0310    By_da0drho_fit = ppval(equil_fit.By_fit.pp_da0drho,rho_fit1);
0311    By_dandrho_fit = ppval(equil_fit.By_fit.pp_dandrho,rho_fit1);
0312    By_dbndrho_fit = ppval(equil_fit.By_fit.pp_dbndrho,rho_fit1);
0313    By_d2a0drho2_fit = ppval(equil_fit.By_fit.pp_d2a0drho2,rho_fit1);
0314    By_d2andrho2_fit = ppval(equil_fit.By_fit.pp_d2andrho2,rho_fit1);
0315    By_d2bndrho2_fit = ppval(equil_fit.By_fit.pp_d2bndrho2,rho_fit1);
0316 %
0317    Bz_a0_fit = ppval(equil_fit.Bz_fit.pp_a0,rho_fit1);
0318    Bz_an_fit = ppval(equil_fit.Bz_fit.pp_an,rho_fit1);
0319    Bz_bn_fit = ppval(equil_fit.Bz_fit.pp_bn,rho_fit1);
0320    Bz_da0drho_fit = ppval(equil_fit.Bz_fit.pp_da0drho,rho_fit1);
0321    Bz_dandrho_fit = ppval(equil_fit.Bz_fit.pp_dandrho,rho_fit1);
0322    Bz_dbndrho_fit = ppval(equil_fit.Bz_fit.pp_dbndrho,rho_fit1);
0323    Bz_d2a0drho2_fit = ppval(equil_fit.Bz_fit.pp_d2a0drho2,rho_fit1);
0324    Bz_d2andrho2_fit = ppval(equil_fit.Bz_fit.pp_d2andrho2,rho_fit1);
0325    Bz_d2bndrho2_fit = ppval(equil_fit.Bz_fit.pp_d2bndrho2,rho_fit1);
0326 %
0327    B_a0_fit = ppval(equil_fit.B_fit.pp_a0,rho_fit1);
0328    B_an_fit = ppval(equil_fit.B_fit.pp_an,rho_fit1);
0329    B_bn_fit = ppval(equil_fit.B_fit.pp_bn,rho_fit1);
0330    B_da0drho_fit = ppval(equil_fit.B_fit.pp_da0drho,rho_fit1);
0331    B_dandrho_fit = ppval(equil_fit.B_fit.pp_dandrho,rho_fit1);
0332    B_dbndrho_fit = ppval(equil_fit.B_fit.pp_dbndrho,rho_fit1);
0333    B_d2a0drho2_fit = ppval(equil_fit.B_fit.pp_d2a0drho2,rho_fit1);
0334    B_d2andrho2_fit = ppval(equil_fit.B_fit.pp_d2andrho2,rho_fit1);
0335    B_d2bndrho2_fit = ppval(equil_fit.B_fit.pp_d2bndrho2,rho_fit1);
0336 %
0337    BP_a0_fit = ppval(equil_fit.BP_fit.pp_a0,rho_fit1);
0338    BP_an_fit = ppval(equil_fit.BP_fit.pp_an,rho_fit1);
0339    BP_bn_fit = ppval(equil_fit.BP_fit.pp_bn,rho_fit1);
0340    BP_da0drho_fit = ppval(equil_fit.BP_fit.pp_da0drho,rho_fit1);
0341    BP_dandrho_fit = ppval(equil_fit.BP_fit.pp_dandrho,rho_fit1);
0342    BP_dbndrho_fit = ppval(equil_fit.BP_fit.pp_dbndrho,rho_fit1);
0343    BP_d2a0drho2_fit = ppval(equil_fit.BP_fit.pp_d2a0drho2,rho_fit1);
0344    BP_d2andrho2_fit = ppval(equil_fit.BP_fit.pp_d2andrho2,rho_fit1);
0345    BP_d2bndrho2_fit = ppval(equil_fit.BP_fit.pp_d2bndrho2,rho_fit1);
0346 %
0347    r_a0_fit = ppval(equil_fit.r_fit.pp_a0,rho_fit1);
0348    r_an_fit = ppval(equil_fit.r_fit.pp_an,rho_fit1);
0349    r_bn_fit = ppval(equil_fit.r_fit.pp_bn,rho_fit1);
0350    r_da0drho_fit = ppval(equil_fit.r_fit.pp_da0drho,rho_fit1);
0351    r_dandrho_fit = ppval(equil_fit.r_fit.pp_dandrho,rho_fit1);
0352    r_dbndrho_fit = ppval(equil_fit.r_fit.pp_dbndrho,rho_fit1); 
0353    r_d2a0drho2_fit = ppval(equil_fit.r_fit.pp_d2a0drho2,rho_fit1);
0354    r_d2andrho2_fit = ppval(equil_fit.r_fit.pp_d2andrho2,rho_fit1);
0355    r_d2bndrho2_fit = ppval(equil_fit.r_fit.pp_d2bndrho2,rho_fit1);   
0356 %
0357    calpha_a0_fit = ppval(equil_fit.calpha_fit.pp_a0,rho_fit1);
0358    calpha_an_fit = ppval(equil_fit.calpha_fit.pp_an,rho_fit1);
0359    calpha_bn_fit = ppval(equil_fit.calpha_fit.pp_bn,rho_fit1);
0360    calpha_da0drho_fit = ppval(equil_fit.calpha_fit.pp_da0drho,rho_fit1);
0361    calpha_dandrho_fit = ppval(equil_fit.calpha_fit.pp_dandrho,rho_fit1);
0362    calpha_dbndrho_fit = ppval(equil_fit.calpha_fit.pp_dbndrho,rho_fit1);
0363    calpha_d2a0drho2_fit = ppval(equil_fit.calpha_fit.pp_d2a0drho2,rho_fit1);
0364    calpha_d2andrho2_fit = ppval(equil_fit.calpha_fit.pp_d2andrho2,rho_fit1);
0365    calpha_d2bndrho2_fit = ppval(equil_fit.calpha_fit.pp_d2bndrho2,rho_fit1);
0366 %
0367    salpha_a0_fit = ppval(equil_fit.salpha_fit.pp_a0,rho_fit1);
0368    salpha_an_fit = ppval(equil_fit.salpha_fit.pp_an,rho_fit1);
0369    salpha_bn_fit = ppval(equil_fit.salpha_fit.pp_bn,rho_fit1);
0370    salpha_da0drho_fit = ppval(equil_fit.salpha_fit.pp_da0drho,rho_fit1);
0371    salpha_dandrho_fit = ppval(equil_fit.salpha_fit.pp_dandrho,rho_fit1);
0372    salpha_dbndrho_fit = ppval(equil_fit.salpha_fit.pp_dbndrho,rho_fit1);
0373    salpha_d2a0drho2_fit = ppval(equil_fit.salpha_fit.pp_d2a0drho2,rho_fit1);
0374    salpha_d2andrho2_fit = ppval(equil_fit.salpha_fit.pp_d2andrho2,rho_fit1);
0375    salpha_d2bndrho2_fit = ppval(equil_fit.salpha_fit.pp_d2bndrho2,rho_fit1);
0376 %
0377    gradrho_a0_fit = ppval(equil_fit.gradrho_fit.pp_a0,rho_fit1);
0378    gradrho_an_fit = ppval(equil_fit.gradrho_fit.pp_an,rho_fit1);
0379    gradrho_bn_fit = ppval(equil_fit.gradrho_fit.pp_bn,rho_fit1);
0380    gradrho_da0drho_fit = ppval(equil_fit.gradrho_fit.pp_da0drho,rho_fit1);
0381    gradrho_dandrho_fit = ppval(equil_fit.gradrho_fit.pp_dandrho,rho_fit1);
0382    gradrho_dbndrho_fit = ppval(equil_fit.gradrho_fit.pp_dbndrho,rho_fit1);
0383    gradrho_d2a0drho2_fit = ppval(equil_fit.gradrho_fit.pp_d2a0drho2,rho_fit1);
0384    gradrho_d2andrho2_fit = ppval(equil_fit.gradrho_fit.pp_d2andrho2,rho_fit1);
0385    gradrho_d2bndrho2_fit = ppval(equil_fit.gradrho_fit.pp_d2bndrho2,rho_fit1);
0386 %
0387 % Calculation of derivatives with home made interpolation using polynomials (spline...) that has a C counterpart
0388 %
0389    x_a0_fit1 = ppval_yp(equil_fit.x_fit.pp_a0,rho_fit1);
0390    x_an_fit1 = ppval_yp(equil_fit.x_fit.pp_an,rho_fit1);
0391    x_bn_fit1 = ppval_yp(equil_fit.x_fit.pp_bn,rho_fit1);
0392    x_da0drho_fit1 = ppval_yp(equil_fit.x_fit.pp_da0drho,rho_fit1);
0393    x_dandrho_fit1 = ppval_yp(equil_fit.x_fit.pp_dandrho,rho_fit1);
0394    x_dbndrho_fit1 = ppval_yp(equil_fit.x_fit.pp_dbndrho,rho_fit1);
0395    x_d2a0drho2_fit1 = ppval_yp(equil_fit.x_fit.pp_d2a0drho2,rho_fit1);
0396    x_d2andrho2_fit1 = ppval_yp(equil_fit.x_fit.pp_d2andrho2,rho_fit1);
0397    x_d2bndrho2_fit1 = ppval_yp(equil_fit.x_fit.pp_d2bndrho2,rho_fit1);
0398 %
0399    y_a0_fit1 = ppval_yp(equil_fit.y_fit.pp_a0,rho_fit1);
0400    y_an_fit1 = ppval_yp(equil_fit.y_fit.pp_an,rho_fit1);
0401    y_bn_fit1 = ppval_yp(equil_fit.y_fit.pp_bn,rho_fit1);
0402    y_da0drho_fit1 = ppval_yp(equil_fit.y_fit.pp_da0drho,rho_fit1);
0403    y_dandrho_fit1 = ppval_yp(equil_fit.y_fit.pp_dandrho,rho_fit1);
0404    y_dbndrho_fit1 = ppval_yp(equil_fit.y_fit.pp_dbndrho,rho_fit1);
0405    y_d2a0drho2_fit1 = ppval_yp(equil_fit.y_fit.pp_d2a0drho2,rho_fit1);
0406    y_d2andrho2_fit1 = ppval_yp(equil_fit.y_fit.pp_d2andrho2,rho_fit1);
0407    y_d2bndrho2_fit1 = ppval_yp(equil_fit.y_fit.pp_d2bndrho2,rho_fit1);
0408 %
0409    Bx_a0_fit1 = ppval_yp(equil_fit.Bx_fit.pp_a0,rho_fit1);
0410    Bx_an_fit1 = ppval_yp(equil_fit.Bx_fit.pp_an,rho_fit1);
0411    Bx_bn_fit1 = ppval_yp(equil_fit.Bx_fit.pp_bn,rho_fit1);
0412    Bx_da0drho_fit1 = ppval_yp(equil_fit.Bx_fit.pp_da0drho,rho_fit1);
0413    Bx_dandrho_fit1 = ppval_yp(equil_fit.Bx_fit.pp_dandrho,rho_fit1);
0414    Bx_dbndrho_fit1 = ppval_yp(equil_fit.Bx_fit.pp_dbndrho,rho_fit1);
0415    Bx_d2a0drho2_fit1 = ppval_yp(equil_fit.Bx_fit.pp_d2a0drho2,rho_fit1);
0416    Bx_d2andrho2_fit1 = ppval_yp(equil_fit.Bx_fit.pp_d2andrho2,rho_fit1);
0417    Bx_d2bndrho2_fit1 = ppval_yp(equil_fit.Bx_fit.pp_d2bndrho2,rho_fit1);
0418 %
0419    By_a0_fit1 = ppval_yp(equil_fit.By_fit.pp_a0,rho_fit1);
0420    By_an_fit1 = ppval_yp(equil_fit.By_fit.pp_an,rho_fit1);
0421    By_bn_fit1 = ppval_yp(equil_fit.By_fit.pp_bn,rho_fit1);
0422    By_da0drho_fit1 = ppval_yp(equil_fit.By_fit.pp_da0drho,rho_fit1);
0423    By_dandrho_fit1 = ppval_yp(equil_fit.By_fit.pp_dandrho,rho_fit1);
0424    By_dbndrho_fit1 = ppval_yp(equil_fit.By_fit.pp_dbndrho,rho_fit1);
0425    By_d2a0drho2_fit1 = ppval_yp(equil_fit.By_fit.pp_d2a0drho2,rho_fit1);
0426    By_d2andrho2_fit1 = ppval_yp(equil_fit.By_fit.pp_d2andrho2,rho_fit1);
0427    By_d2bndrho2_fit1 = ppval_yp(equil_fit.By_fit.pp_d2bndrho2,rho_fit1);
0428 %
0429    Bz_a0_fit1 = ppval_yp(equil_fit.Bz_fit.pp_a0,rho_fit1);
0430    Bz_an_fit1 = ppval_yp(equil_fit.Bz_fit.pp_an,rho_fit1);
0431    Bz_bn_fit1 = ppval_yp(equil_fit.Bz_fit.pp_bn,rho_fit1);
0432    Bz_da0drho_fit1 = ppval_yp(equil_fit.Bz_fit.pp_da0drho,rho_fit1);
0433    Bz_dandrho_fit1 = ppval_yp(equil_fit.Bz_fit.pp_dandrho,rho_fit1);
0434    Bz_dbndrho_fit1 = ppval_yp(equil_fit.Bz_fit.pp_dbndrho,rho_fit1);
0435    Bz_d2a0drho2_fit1 = ppval_yp(equil_fit.Bz_fit.pp_d2a0drho2,rho_fit1);
0436    Bz_d2andrho2_fit1 = ppval_yp(equil_fit.Bz_fit.pp_d2andrho2,rho_fit1);
0437    Bz_d2bndrho2_fit1 = ppval_yp(equil_fit.Bz_fit.pp_d2bndrho2,rho_fit1);
0438 %
0439    B_a0_fit1 = ppval_yp(equil_fit.B_fit.pp_a0,rho_fit1);
0440    B_an_fit1 = ppval_yp(equil_fit.B_fit.pp_an,rho_fit1);
0441    B_bn_fit1 = ppval_yp(equil_fit.B_fit.pp_bn,rho_fit1);
0442    B_da0drho_fit1 = ppval_yp(equil_fit.B_fit.pp_da0drho,rho_fit1);
0443    B_dandrho_fit1 = ppval_yp(equil_fit.B_fit.pp_dandrho,rho_fit1);
0444    B_dbndrho_fit1 = ppval_yp(equil_fit.B_fit.pp_dbndrho,rho_fit1);
0445    B_d2a0drho2_fit1 = ppval_yp(equil_fit.B_fit.pp_d2a0drho2,rho_fit1);
0446    B_d2andrho2_fit1 = ppval_yp(equil_fit.B_fit.pp_d2andrho2,rho_fit1);
0447    B_d2bndrho2_fit1 = ppval_yp(equil_fit.B_fit.pp_d2bndrho2,rho_fit1);
0448 %
0449    BP_a0_fit1 = ppval_yp(equil_fit.BP_fit.pp_a0,rho_fit1);
0450    BP_an_fit1 = ppval_yp(equil_fit.BP_fit.pp_an,rho_fit1);
0451    BP_bn_fit1 = ppval_yp(equil_fit.BP_fit.pp_bn,rho_fit1);
0452    BP_da0drho_fit1 = ppval_yp(equil_fit.BP_fit.pp_da0drho,rho_fit1);
0453    BP_dandrho_fit1 = ppval_yp(equil_fit.BP_fit.pp_dandrho,rho_fit1);
0454    BP_dbndrho_fit1 = ppval_yp(equil_fit.BP_fit.pp_dbndrho,rho_fit1);
0455    BP_d2a0drho2_fit1 = ppval_yp(equil_fit.BP_fit.pp_d2a0drho2,rho_fit1);
0456    BP_d2andrho2_fit1 = ppval_yp(equil_fit.BP_fit.pp_d2andrho2,rho_fit1);
0457    BP_d2bndrho2_fit1 = ppval_yp(equil_fit.BP_fit.pp_d2bndrho2,rho_fit1);
0458 %
0459    r_a0_fit1 = ppval_yp(equil_fit.r_fit.pp_a0,rho_fit1);
0460    r_an_fit1 = ppval_yp(equil_fit.r_fit.pp_an,rho_fit1);
0461    r_bn_fit1 = ppval_yp(equil_fit.r_fit.pp_bn,rho_fit1);
0462    r_da0drho_fit1 = ppval_yp(equil_fit.r_fit.pp_da0drho,rho_fit1);
0463    r_dandrho_fit1 = ppval_yp(equil_fit.r_fit.pp_dandrho,rho_fit1);
0464    r_dbndrho_fit1 = ppval_yp(equil_fit.r_fit.pp_dbndrho,rho_fit1);
0465    r_d2a0drho2_fit1 = ppval_yp(equil_fit.r_fit.pp_d2a0drho2,rho_fit1);
0466    r_d2andrho2_fit1 = ppval_yp(equil_fit.r_fit.pp_d2andrho2,rho_fit1);
0467    r_d2bndrho2_fit1 = ppval_yp(equil_fit.r_fit.pp_d2bndrho2,rho_fit1);
0468 %
0469    calpha_a0_fit1 = ppval_yp(equil_fit.calpha_fit.pp_a0,rho_fit1);
0470    calpha_an_fit1 = ppval_yp(equil_fit.calpha_fit.pp_an,rho_fit1);
0471    calpha_bn_fit1 = ppval_yp(equil_fit.calpha_fit.pp_bn,rho_fit1);
0472    calpha_da0drho_fit1 = ppval_yp(equil_fit.calpha_fit.pp_da0drho,rho_fit1);
0473    calpha_dandrho_fit1 = ppval_yp(equil_fit.calpha_fit.pp_dandrho,rho_fit1);
0474    calpha_dbndrho_fit1 = ppval_yp(equil_fit.calpha_fit.pp_dbndrho,rho_fit1);
0475    calpha_d2a0drho2_fit1 = ppval_yp(equil_fit.calpha_fit.pp_d2a0drho2,rho_fit1);
0476    calpha_d2andrho2_fit1 = ppval_yp(equil_fit.calpha_fit.pp_d2andrho2,rho_fit1);
0477    calpha_d2bndrho2_fit1 = ppval_yp(equil_fit.calpha_fit.pp_d2bndrho2,rho_fit1);
0478 %
0479    salpha_a0_fit1 = ppval_yp(equil_fit.salpha_fit.pp_a0,rho_fit1);
0480    salpha_an_fit1 = ppval_yp(equil_fit.salpha_fit.pp_an,rho_fit1);
0481    salpha_bn_fit1 = ppval_yp(equil_fit.salpha_fit.pp_bn,rho_fit1);
0482    salpha_da0drho_fit1 = ppval_yp(equil_fit.salpha_fit.pp_da0drho,rho_fit1);
0483    salpha_dandrho_fit1 = ppval_yp(equil_fit.salpha_fit.pp_dandrho,rho_fit1);
0484    salpha_dbndrho_fit1 = ppval_yp(equil_fit.salpha_fit.pp_dbndrho,rho_fit1);
0485    salpha_d2a0drho2_fit1 = ppval_yp(equil_fit.salpha_fit.pp_d2a0drho2,rho_fit1);
0486    salpha_d2andrho2_fit1 = ppval_yp(equil_fit.salpha_fit.pp_d2andrho2,rho_fit1);
0487    salpha_d2bndrho2_fit1 = ppval_yp(equil_fit.salpha_fit.pp_d2bndrho2,rho_fit1);
0488 %
0489    gradrho_a0_fit1 = ppval_yp(equil_fit.gradrho_fit.pp_a0,rho_fit1);
0490    gradrho_an_fit1 = ppval_yp(equil_fit.gradrho_fit.pp_an,rho_fit1);
0491    gradrho_bn_fit1 = ppval_yp(equil_fit.gradrho_fit.pp_bn,rho_fit1);
0492    gradrho_da0drho_fit1 = ppval_yp(equil_fit.gradrho_fit.pp_da0drho,rho_fit1);
0493    gradrho_dandrho_fit1 = ppval_yp(equil_fit.gradrho_fit.pp_dandrho,rho_fit1);
0494    gradrho_dbndrho_fit1 = ppval_yp(equil_fit.gradrho_fit.pp_dbndrho,rho_fit1);
0495    gradrho_d2a0drho2_fit1 = ppval_yp(equil_fit.gradrho_fit.pp_d2a0drho2,rho_fit1);
0496    gradrho_d2andrho2_fit1 = ppval_yp(equil_fit.gradrho_fit.pp_d2andrho2,rho_fit1);
0497    gradrho_d2bndrho2_fit1 = ppval_yp(equil_fit.gradrho_fit.pp_d2bndrho2,rho_fit1);
0498 %
0499 % Build interpolated quantities
0500 %
0501    for ii = 1:ntheta_fit, 
0502       [x_fit(:,ii),dxdtheta_fit(:,ii),d2xdtheta2_fit(:,ii),dxdrho_fit(:,ii),d2xdthetadrho_fit(:,ii),d2xdrho2_fit(:,ii)] = calcval_yp(equil_fit,theta_fit(ii),x_a0_fit(:),x_an_fit,x_bn_fit,x_da0drho_fit(:),x_dandrho_fit,x_dbndrho_fit,x_d2a0drho2_fit(:),x_d2andrho2_fit,x_d2bndrho2_fit);
0503       [y_fit(:,ii),dydtheta_fit(:,ii),d2ydtheta2_fit(:,ii),dydrho_fit(:,ii),d2ydthetadrho_fit(:,ii),d2ydrho2_fit(:,ii)] = calcval_yp(equil_fit,theta_fit(ii),y_a0_fit(:),y_an_fit,y_bn_fit,y_da0drho_fit(:),y_dandrho_fit,y_dbndrho_fit,y_d2a0drho2_fit(:),y_d2andrho2_fit,y_d2bndrho2_fit);
0504       %
0505       [r_fit(:,ii),drdtheta_fit(:,ii),d2rdtheta2_fit(:,ii),drdrho_fit(:,ii),d2rdthetadrho_fit(:,ii),d2rdrho2_fit(:,ii)] = calcval_yp(equil_fit,theta_fit(ii),r_a0_fit(:),r_an_fit,r_bn_fit,r_da0drho_fit(:),r_dandrho_fit,r_dbndrho_fit,r_d2a0drho2_fit(:),r_d2andrho2_fit,r_d2bndrho2_fit);
0506       %
0507       [Bx_fit(:,ii),dBxdtheta_fit(:,ii),d2Bxdtheta2_fit(:,ii),dBxdrho_fit(:,ii),d2Bxdthetadrho_fit(:,ii),d2Bxdrho2_fit(:,ii)] = calcval_yp(equil_fit,theta_fit(ii),Bx_a0_fit(:),Bx_an_fit,Bx_bn_fit,Bx_da0drho_fit(:),Bx_dandrho_fit,Bx_dbndrho_fit,Bx_d2a0drho2_fit(:),Bx_d2andrho2_fit,Bx_d2bndrho2_fit);
0508       [By_fit(:,ii),dBydtheta_fit(:,ii),d2Bydtheta2_fit(:,ii),dBydrho_fit(:,ii),d2Bydthetadrho_fit(:,ii),d2Bydrho2_fit(:,ii)] = calcval_yp(equil_fit,theta_fit(ii),By_a0_fit(:),By_an_fit,By_bn_fit,By_da0drho_fit(:),By_dandrho_fit,By_dbndrho_fit,By_d2a0drho2_fit(:),By_d2andrho2_fit,By_d2bndrho2_fit);
0509       [Bz_fit(:,ii),dBzdtheta_fit(:,ii),d2Bzdtheta2_fit(:,ii),dBzdrho_fit(:,ii),d2Bzdthetadrho_fit(:,ii),d2Bzdrho2_fit(:,ii)] = calcval_yp(equil_fit,theta_fit(ii),Bz_a0_fit(:),Bz_an_fit,Bz_bn_fit,Bz_da0drho_fit(:),Bz_dandrho_fit,Bz_dbndrho_fit,Bz_d2a0drho2_fit(:),Bz_d2andrho2_fit,Bz_d2bndrho2_fit);
0510       %
0511       [B_fit(:,ii),dBdtheta_fit(:,ii),d2Bdtheta2_fit(:,ii),dBdrho_fit(:,ii),d2Bdthetadrho_fit(:,ii),d2Bdrho2_fit(:,ii)] = calcval_yp(equil_fit,theta_fit(ii),B_a0_fit(:),B_an_fit,B_bn_fit,B_da0drho_fit(:),B_dandrho_fit,B_dbndrho_fit,B_d2a0drho2_fit(:),B_d2andrho2_fit,B_d2bndrho2_fit);
0512       [BP_fit(:,ii),dBPdtheta_fit(:,ii),d2BPdtheta2_fit(:,ii),dBPdrho_fit(:,ii),d2BPdthetadrho_fit(:,ii),d2BPdrho2_fit(:,ii)] = calcval_yp(equil_fit,theta_fit(ii),BP_a0_fit(:),BP_an_fit,BP_bn_fit,BP_da0drho_fit(:),BP_dandrho_fit,BP_dbndrho_fit,BP_d2a0drho2_fit(:),BP_d2andrho2_fit,BP_d2bndrho2_fit);
0513       %
0514       [calpha_fit(:,ii),dcalphadtheta_fit(:,ii),d2calphadtheta2_fit(:,ii),dcalphadrho_fit(:,ii),d2calphadthetadrho_fit(:,ii),d2calphadrho2_fit(:,ii)] = calcval_yp(equil_fit,theta_fit(ii),calpha_a0_fit(:),calpha_an_fit,calpha_bn_fit,calpha_da0drho_fit(:),calpha_dandrho_fit,calpha_dbndrho_fit,calpha_d2a0drho2_fit(:),calpha_d2andrho2_fit,calpha_d2bndrho2_fit);
0515       [salpha_fit(:,ii),dsalphadtheta_fit(:,ii),d2salphadtheta2_fit(:,ii),dsalphadrho_fit(:,ii),d2salphadthetadrho_fit(:,ii),d2salphadrho2_fit(:,ii)] = calcval_yp(equil_fit,theta_fit(ii),salpha_a0_fit(:),salpha_an_fit,salpha_bn_fit,salpha_da0drho_fit(:),salpha_dandrho_fit,salpha_dbndrho_fit,salpha_d2a0drho2_fit(:),salpha_d2andrho2_fit,salpha_d2bndrho2_fit);     
0516       %
0517       [gradrho_fit(:,ii),dgradrhodtheta_fit(:,ii),d2gradrhodtheta2_fit(:,ii),dgradrhodrho_fit(:,ii),d2gradrhodthetadrho_fit(:,ii),d2gradrhodrho2_fit(:,ii)] = calcval_yp(equil_fit,theta_fit(ii),gradrho_a0_fit(:),gradrho_an_fit,gradrho_bn_fit,gradrho_da0drho_fit(:),gradrho_dandrho_fit,gradrho_dbndrho_fit,gradrho_d2a0drho2_fit(:),gradrho_d2andrho2_fit,gradrho_d2bndrho2_fit);     
0518       %
0519       if isinf(Rp),
0520          Upsilon_fit = ones(size(x_fit(:,ii)));
0521          dUpsilondrho_fit = zeros(size(x_fit(:,ii)));
0522          dUpsilondtheta_fit = zeros(size(x_fit(:,ii)));
0523          d2Upsilondtheta2_fit = zeros(size(x_fit(:,ii)));
0524          d2Upsilondthetadrho_fit = zeros(size(x_fit(:,ii)));
0525          d2Upsilondrho2_fit = zeros(size(x_fit(:,ii)));
0526       else
0527          Upsilon_fit = 1 + x_fit(:,ii)/Rp;
0528          dUpsilondrho_fit = dxdrho_fit(:,ii)/Rp;
0529          dUpsilondtheta_fit =  dxdtheta_fit(:,ii)/Rp;
0530          d2Upsilondtheta2_fit = d2xdtheta2_fit(:,ii)/Rp;
0531          d2Upsilondthetadrho_fit = d2xdthetadrho_fit(:,ii)/Rp;
0532          d2Upsilondrho2_fit = d2xdrho2_fit(:,ii)/Rp;
0533 
0534       end
0535       %
0536       [x_fit1(:,ii),dxdtheta_fit1(:,ii),d2xdtheta2_fit1(:,ii),dxdrho_fit1(:,ii),d2xdthetadrho_fit1(:,ii),d2xdrho2_fit1(:,ii)] = calcval_yp(equil_fit,theta_fit(ii),x_a0_fit1,x_an_fit1,x_bn_fit,x_da0drho_fit1,x_dandrho_fit1,x_dbndrho_fit1,x_d2a0drho2_fit1,x_d2andrho2_fit1,x_d2bndrho2_fit1);
0537       [y_fit1(:,ii),dydtheta_fit1(:,ii),d2ydtheta2_fit1(:,ii),dydrho_fit1(:,ii),d2ydthetadrho_fit1(:,ii),d2ydrho2_fit1(:,ii)] = calcval_yp(equil_fit,theta_fit(ii),y_a0_fit1,y_an_fit1,y_bn_fit,y_da0drho_fit1,y_dandrho_fit1,y_dbndrho_fit1,y_d2a0drho2_fit1,y_d2andrho2_fit1,y_d2bndrho2_fit1);
0538       %
0539       [r_fit1(:,ii),drdtheta_fit1(:,ii),d2rdtheta2_fit1(:,ii),drdrho_fit1(:,ii),d2rdthetadrho_fit1(:,ii),d2rdrho2_fit1(:,ii)] = calcval_yp(equil_fit,theta_fit(ii),r_a0_fit1,r_an_fit1,r_bn_fit,r_da0drho_fit1,r_dandrho_fit1,r_dbndrho_fit1,r_d2a0drho2_fit1,r_d2andrho2_fit1,r_d2bndrho2_fit1);      
0540       %
0541       [Bx_fit1(:,ii),dBxdtheta_fit1(:,ii),d2Bxdtheta2_fit1(:,ii),dBxdrho_fit1(:,ii),d2Bxdthetadrho_fit1(:,ii),d2Bxdrho2_fit1(:,ii)] = calcval_yp(equil_fit,theta_fit(ii),Bx_a0_fit1,Bx_an_fit1,Bx_bn_fit1,Bx_da0drho_fit1,Bx_dandrho_fit1,Bx_dbndrho_fit1,Bx_d2a0drho2_fit1,Bx_d2andrho2_fit1,Bx_d2bndrho2_fit1);
0542       [By_fit1(:,ii),dBydtheta_fit1(:,ii),d2Bydtheta2_fit1(:,ii),dBydrho_fit1(:,ii),d2Bydthetadrho_fit1(:,ii),d2Bxydrho2_fit1(:,ii)] = calcval_yp(equil_fit,theta_fit(ii),By_a0_fit1,By_an_fit1,By_bn_fit1,By_da0drho_fit1,By_dandrho_fit1,By_dbndrho_fit1,By_d2a0drho2_fit1,By_d2andrho2_fit1,By_d2bndrho2_fit1);
0543       [Bz_fit1(:,ii),dBzdtheta_fit1(:,ii),d2Bzdtheta2_fit1(:,ii),dBzdrho_fit1(:,ii),d2Bzdthetadrho_fit1(:,ii),d2Bzdrho2_fit1(:,ii)] = calcval_yp(equil_fit,theta_fit(ii),Bz_a0_fit1,Bz_an_fit1,Bz_bn_fit1,Bz_da0drho_fit1,Bz_dandrho_fit1,Bz_dbndrho_fit1,Bz_d2a0drho2_fit1,Bz_d2andrho2_fit1,Bz_d2bndrho2_fit1);
0544       %
0545       [B_fit1(:,ii),dBdtheta_fit1(:,ii),d2Bdtheta2_fit1(:,ii),dBdrho_fit1(:,ii),d2Bdthetadrho_fit1(:,ii),d2Bdrho2_fit1(:,ii)] = calcval_yp(equil_fit,theta_fit(ii),B_a0_fit1,B_an_fit1,B_bn_fit1,B_da0drho_fit1,B_dandrho_fit1,B_dbndrho_fit1,B_d2a0drho2_fit1,B_d2andrho2_fit1,B_d2bndrho2_fit1);
0546       [BP_fit1(:,ii),dBPdtheta_fit1(:,ii),d2BPdtheta2_fit1(:,ii),dBPdrho_fit1(:,ii),d2BPdthetadrho_fit1(:,ii),d2BPdrho2_fit1(:,ii)] = calcval_yp(equil_fit,theta_fit(ii),BP_a0_fit1,BP_an_fit1,BP_bn_fit1,BP_da0drho_fit1,BP_dandrho_fit1,BP_dbndrho_fit1,BP_d2a0drho2_fit1,BP_d2andrho2_fit1,BP_d2bndrho2_fit1);
0547       %
0548       [calpha_fit1(:,ii),dcalphadtheta_fit1(:,ii),d2calphadtheta2_fit1(:,ii),dcalphadrho_fit1(:,ii),d2calphadthetadrho_fit1(:,ii),d2calphadrho2_fit1(:,ii)] = calcval_yp(equil_fit,theta_fit(ii),calpha_a0_fit1,calpha_an_fit1,calpha_bn_fit1,calpha_da0drho_fit1,calpha_dandrho_fit1,calpha_dbndrho_fit1,calpha_d2a0drho2_fit1,calpha_d2andrho2_fit1,calpha_d2bndrho2_fit1);
0549       [salpha_fit1(:,ii),dsalphadtheta_fit1(:,ii),d2salphadtheta2_fit1(:,ii),dsalphadrho_fit1(:,ii),d2salphadthetadrho_fit1(:,ii),d2salphadrho2_fit1(:,ii)] = calcval_yp(equil_fit,theta_fit(ii),salpha_a0_fit1,salpha_an_fit1,salpha_bn_fit1,salpha_da0drho_fit1,salpha_dandrho_fit1,salpha_dbndrho_fit1,salpha_d2a0drho2_fit1,salpha_d2andrho2_fit1,salpha_d2bndrho2_fit1);
0550       %
0551       [gradrho_fit1(:,ii),dgradrhodtheta_fit1(:,ii),d2gradrhodtheta2_fit1(:,ii),dgradrhodrho_fit1(:,ii),d2gradrhodthetadrho_fit1(:,ii),d2gradrhodrho2_fit1(:,ii)] = calcval_yp(equil_fit,theta_fit(ii),gradrho_a0_fit1(:),gradrho_an_fit1,gradrho_bn_fit1,gradrho_da0drho_fit1(:),gradrho_dandrho_fit1,gradrho_dbndrho_fit1,gradrho_d2a0drho2_fit1,gradrho_d2andrho2_fit1,gradrho_d2bndrho2_fit1);     
0552       %
0553       if isinf(Rp),
0554          Upsilon_fit1 = ones(size(x_fit1(:,ii)));
0555          dUpsilondrho_fit1 = zeros(size(x_fit1(:,ii)));
0556          dUpsilondtheta_fit1 = zeros(size(x_fit1(:,ii)));
0557          d2Upsilondtheta2_fit1 = zeros(size(x_fit1(:,ii)));
0558          d2Upsilondthetadrho_fit1 = zeros(size(x_fit1(:,ii)));
0559          d2Upsilondrho2_fit1 = zeros(size(x_fit1(:,ii)));
0560       else
0561          Upsilon_fit1 = 1 + x_fit1(:,ii)/Rp;
0562          dUpsilondrho_fit1 = dxdrho_fit1(:,ii)/Rp;
0563          dUpsilondtheta_fit1 =  dxdtheta_fit1(:,ii)/Rp;
0564          d2Upsilondtheta2_fit1 = d2xdtheta2_fit1(:,ii)/Rp;
0565          d2Upsilondthetadrho_fit1 = d2xdthetadrho_fit1(:,ii)/Rp;
0566          d2Upsilondrho2_fit1 = d2xdrho2_fit1(:,ii)/Rp;
0567       end
0568       %
0569    end
0570    %
0571    % Display results
0572    %
0573    scale_factor = 1.3;
0574    figure('Name','Magnetic equilibrium profile parameters'),clf
0575    set(gcf,'units','inches','position',[0.5,1,10*scale_factor,6*scale_factor]);
0576    set(gca,'position',[0,0,1,1]);
0577    %
0578    subplot('position',[1/12,0.5+1/16,5/24,5/16])
0579    [ax] = graph1D_jd(rho_fit,psin_fit,0,0,'\rho','\psi_n','',NaN,[0,1],[-eps,1],'-','none','b',2);%ppval matlab
0580    [ax] = graph1D_jd(rho_fit,psin_fit1,0,0,'\rho','\psi_n','',NaN,[0,1],[-eps,1],'-','none','c',1);%ppval_yp matlab
0581    [ax] = graph1D_jd(prho(1:nr:end),psin(1:nr:end),0,0,'\rho','\psi_n','',NaN,[0,1],[-eps,1],'--','non','r',1);%numerical interpolation
0582    if isfield(equil_fit,'Te_fit'),
0583       subplot('position',[1/3+1/12,0.5+1/16,5/24,5/16])
0584       [ax] = graph1D_jd(rho_fit,Te_fit,0,0,'\psi','T_e (keV)',strrep(equil_in.id,'_','-'),NaN,[0,1],[-eps,max(equil_in.pTe)],'-','none','b',2);
0585       [ax] = graph1D_jd(rho_fit,Te_fit1,0,0,'\psi','T_e (keV)',strrep(equil_in.id,'_','-'),NaN,[0,1],[-eps,max(equil_in.pTe)],'-','none','c',1);
0586       [ax] = graph1D_jd(prho(1:nr:end),pTe(1:nr:end),0,0,'\rho','T_e (keV)',[strrep(equil_in.id,'_','-'),', blue: ppval, cyan: ppval-yp, red: numerical'],NaN,[0,1],[-eps,max(equil_in.pTe)],'--','non','r',1);
0587    end
0588    if isfield(equil_fit,'ne_fit'),
0589       subplot('position',[2/3+1/12,0.5+1/16,5/24,5/16])
0590       [ax] = graph1D_jd(rho_fit,ne_fit/1e19,0,0,'\rho','n_e (10^{19} m^{-3})','',NaN,[0,1],[-eps,max(equil_in.pne/1e19)],'-','none','b',2);
0591       [ax] = graph1D_jd(rho_fit,ne_fit1/1e19,0,0,'\rho','n_e (10^{19} m^{-3})','',NaN,[0,1],[-eps,max(equil_in.pne/1e19)],'-','none','c',1);
0592       [ax] = graph1D_jd(prho(1:nr:end),pne(1:nr:end)/1e19,0,0,'\rho','n_e (10^{19} m^{-3})','',NaN,[0,1],[-eps,max(equil_in.pne/1e19)],'--','non','r',1);
0593    end
0594    if isfield(equil_fit,'zTi_fit'),
0595       subplot('position',[1/12,1/8,5/24,5/16])
0596       [ax] = graph1D_jd(rho_fit,zTi_fit,0,0,'\rho','T_i (keV)','',NaN,[0,1],[-eps,max(max(equil_in.pzTi))],'-','none','b',2);
0597       [ax] = graph1D_jd(rho_fit,zTi_fit1,0,0,'\rho','T_i (keV)','',NaN,[0,1],[-eps,max(max(equil_in.pzTi))],'-','none','c',1);
0598       [ax] = graph1D_jd(prho(1:nr:end),pzTi(1:nr:end,:),0,0,'\rho','T_i (keV)','',NaN,[0,1],[-eps,max(max(equil_in.pzTi))],'--','non','r',1);
0599    end
0600    if isfield(equil_fit,'zni_fit'),
0601       subplot('position',[1/3+1/12,1/8,5/24,5/16])
0602       [ax] = graph1D_jd(rho_fit,zni_fit/1e19,0,0,'\rho','n_i (10^{19} m^{-3})','',NaN,[0,1],[-eps,max(max(equil_in.pzni/1e19))],'-','none','b',2);
0603       [ax] = graph1D_jd(rho_fit,zni_fit1/1e19,0,0,'\rho','n_i (10^{19} m^{-3})','',NaN,[0,1],[-eps,max(max(equil_in.pzni/1e19))],'-','none','c',1);
0604       [ax] = graph1D_jd(prho(1:nr:end),pzni(1:nr:end,:)/1e19,0,0,'\rho','n_i (10^{19} m^{-3})','',NaN,[0,1],[-eps,max(max(equil_in.pzni/1e19))],'--','non','r',1);
0605       subplot('position',[2/3+1/12,1/8,5/24,5/16])
0606    end
0607    %
0608    figure('Name','1st order derivatives of the magnetic equilibrium profile parameters'),clf
0609    set(gcf,'units','inches','position',[0.5,1,10*scale_factor,6*scale_factor]);
0610    set(gca,'position',[0,0,1,1]);
0611    %
0612    subplot('position',[1/12,0.5+1/16,5/24,5/16])
0613    [ax] = graph1D_jd(rho_fit,dpsindrho_fit,0,0,'','','',NaN,[0,1],[-eps,1],'-','none','b',2);
0614    [ax] = graph1D_jd(rho_fit,dpsindrho_fit1,0,0,'','','',NaN,[0,1],[-eps,1],'-','none','c',1);
0615    [ax] = graph1D_jd(prho(1:nr:end),dpsindrho(1:nr:end),0,0,'\rho','d\psi_nd/\rho','',NaN,[0,1],[min(dpsindrho_fit1),max(dpsindrho_fit1)],'--','non','r',1);
0616    if isfield(equil_fit,'Te_fit'),
0617       subplot('position',[1/3+1/12,0.5+1/16,5/24,5/16])
0618       [ax] = graph1D_jd(rho_fit,dTedrho_fit,0,0,'','',strrep(equil_in.id,'_','-'),NaN,NaN,NaN,'-','none','b',2);
0619       [ax] = graph1D_jd(rho_fit,dTedrho_fit1,0,0,'','',[strrep(equil_in.id,'_','-'),', blue: ppval, cyan: ppval-yp, red: numerical'],NaN,NaN,NaN,'-','none','c',1);
0620       [ax] = graph1D_jd(prho(1:nr:end),dTedrho(1:nr:end),0,0,'\rho','dT_ed/\rho (keV)',[strrep(equil_in.id,'_','-'),', blue: ppval, cyan: ppval-yp, red: numerical'],NaN,NaN,NaN,'--','non','r',1);
0621    end
0622    if isfield(equil_fit,'ne_fit'),
0623       subplot('position',[2/3+1/12,0.5+1/16,5/24,5/16])
0624       [ax] = graph1D_jd(rho_fit,dnedrho_fit/1e19,0,0,'','','',NaN,NaN,NaN,'-','none','b',2);
0625       [ax] = graph1D_jd(rho_fit,dnedrho_fit1/1e19,0,0,'','','',NaN,NaN,NaN,'-','none','c',1);
0626       [ax] = graph1D_jd(prho(1:nr:end),dnedrho(1:nr:end)/1e19,0,0,'\rho','dn_e/d\rho (10^{19} m^{-3})','',NaN,NaN,NaN,'--','non','r',1);
0627    end
0628    if isfield(equil_fit,'zTi_fit'),
0629       subplot('position',[1/12,1/8,5/24,5/16])
0630       [ax] = graph1D_jd(rho_fit,zdTidrho_fit,0,0,'','','',NaN,NaN,NaN,'-','none','b',2);
0631       [ax] = graph1D_jd(rho_fit,zdTidrho_fit1,0,0,'','','',NaN,NaN,NaN,'-','none','c',1);
0632       [ax] = graph1D_jd(prho(1:nr:end),zdTidrho(1:nr:end,:),0,0,'\rho','dT_i/d\rho (keV)','',NaN,NaN,NaN,'--','non','r',1);
0633    end
0634    if isfield(equil_fit,'zni_fit'),
0635       subplot('position',[1/3+1/12,1/8,5/24,5/16])
0636       [ax] = graph1D_jd(rho_fit,zdnidrho_fit/1e19,0,0,'','','',NaN,NaN,NaN,'-','none','b',2);
0637       [ax] = graph1D_jd(rho_fit,zdnidrho_fit1/1e19,0,0,'','','',NaN,NaN,NaN,'-','none','c',1);
0638       [ax] = graph1D_jd(prho(1:nr:end),zdnidrho(1:nr:end,:)/1e19,0,0,'\rho','dn_i/d\rho (10^{19} m^{-3})','',NaN,NaN,NaN,'--','non','r',1);
0639       subplot('position',[2/3+1/12,1/8,5/24,5/16])
0640    end
0641    %
0642    figure('Name','2nd order derivatives of the magnetic equilibrium profile parameters'),clf
0643    set(gcf,'units','inches','position',[0.5,1,10*scale_factor,6*scale_factor]);
0644    set(gca,'position',[0,0,1,1]);
0645    %
0646    subplot('position',[1/12,0.5+1/16,5/24,5/16])
0647    [ax] = graph1D_jd(rho_fit,d2psindrho2_fit,0,0,'','','',NaN,[0,1],[-eps,1],'-','none','b',2);
0648    [ax] = graph1D_jd(rho_fit,d2psindrho2_fit1,0,0,'','','',NaN,[0,1],[-eps,1],'-','none','c',1);
0649    [ax] = graph1D_jd(prho(1:nr:end),d2psindrho2(1:nr:end),0,0,'\rho','d2\psi_nd/\rho2','',NaN,[0,1],[min(d2psindrho2_fit1),max(d2psindrho2_fit1)],'--','non','r',1);
0650    if isfield(equil_fit,'Te_fit'),
0651       subplot('position',[1/3+1/12,0.5+1/16,5/24,5/16])
0652       [ax] = graph1D_jd(rho_fit,d2Tedrho2_fit,0,0,'','',strrep(equil_in.id,'_','-'),NaN,NaN,NaN,'-','none','b',2);
0653       [ax] = graph1D_jd(rho_fit,d2Tedrho2_fit1,0,0,'','',[strrep(equil_in.id,'_','-'),', blue: ppval, cyan: ppval-yp, red: numerical'],NaN,NaN,NaN,'-','none','c',1);
0654       [ax] = graph1D_jd(prho(1:nr:end),d2Tedrho2(1:nr:end),0,0,'\rho','d2T_ed/\rho2 (keV)',[strrep(equil_in.id,'_','-'),', blue: ppval, cyan: ppval-yp, red: numerical'],NaN,NaN,NaN,'--','non','r',1);
0655    end
0656    if isfield(equil_fit,'ne_fit'),
0657       subplot('position',[2/3+1/12,0.5+1/16,5/24,5/16])
0658       [ax] = graph1D_jd(rho_fit,d2nedrho2_fit/1e19,0,0,'','','',NaN,NaN,NaN,'-','none','b',2);
0659       [ax] = graph1D_jd(rho_fit,d2nedrho2_fit1/1e19,0,0,'','','',NaN,NaN,NaN,'-','none','c',1);
0660       [ax] = graph1D_jd(prho(1:nr:end),d2nedrho2(1:nr:end)/1e19,0,0,'\rho','d2n_e/d\rho2 (10^{19} m^{-3})','',NaN,NaN,NaN,'--','non','r',1);
0661    end
0662    if isfield(equil_fit,'zTi_fit'),
0663       subplot('position',[1/12,1/8,5/24,5/16])
0664       [ax] = graph1D_jd(rho_fit,zd2Tidrho2_fit,0,0,'','','',NaN,NaN,NaN,'-','none','b',2);
0665       [ax] = graph1D_jd(rho_fit,zd2Tidrho2_fit1,0,0,'','','',NaN,NaN,NaN,'-','none','c',1);
0666       [ax] = graph1D_jd(prho(1:nr:end),zd2Tidrho2(1:nr:end,:),0,0,'\rho','d2T_i/d\rho2 (keV)','',NaN,NaN,NaN,'--','non','r',1);
0667    end
0668    if isfield(equil_fit,'zni_fit'),
0669       subplot('position',[1/3+1/12,1/8,5/24,5/16])
0670       [ax] = graph1D_jd(rho_fit,zd2nidrho2_fit/1e19,0,0,'','','',NaN,NaN,NaN,'-','none','b',2);
0671       [ax] = graph1D_jd(rho_fit,zd2nidrho2_fit1/1e19,0,0,'','','',NaN,NaN,NaN,'-','none','c',1);
0672       [ax] = graph1D_jd(prho(1:nr:end),zd2nidrho2(1:nr:end,:)/1e19,0,0,'','d2n_i/d\rho2 (10^{19} m^{-3})','',NaN,NaN,NaN,'--','non','r',1);
0673       subplot('position',[2/3+1/12,1/8,5/24,5/16])
0674    end
0675    %
0676    % Magnetic equilibrium coordinates
0677    %
0678    figure('Name','x (majord radius) coordinate'),
0679    [ax] = graph1D_jd(theta_fit,x_fit,0,0,'','','',NaN,NaN,NaN,'-','none','b',2);
0680    [ax] = graph1D_jd(theta_fit,x_fit1,0,0,'','','',NaN,NaN,NaN,'-','none','c',1);
0681    [ax] = graph1D_jd(equil_in.theta(1:2:end),x(1:nr:end,1:2:end),0,0,'\theta (rd)','x (m)',[strrep(equil_in.id,'_','-'),', blue: ppval, cyan: ppval-yp, red: numerical'],NaN,NaN,NaN,'--','non','r',1);
0682    %
0683    figure('Name','dx/dtheta'),
0684    [ax] = graph1D_jd(theta_fit,dxdtheta_fit,0,0,'','','',NaN,NaN,NaN,'-','none','b',2);
0685    [ax] = graph1D_jd(theta_fit,dxdtheta_fit1,0,0,'','','',NaN,NaN,NaN,'-','none','c',1);
0686    [ax] = graph1D_jd(equil_in.theta(1:2:end),dxdtheta(1:nr:end,1:2:end),0,0,'\theta (rd)','dx/d\theta',[strrep(equil_in.id,'_','-'),', blue: ppval, cyan: ppval-yp, red: numerical'],NaN,NaN,NaN,'--','non','r',1);
0687    %
0688    figure('Name','dx/drho'),
0689    [ax] = graph1D_jd(theta_fit,dxdrho_fit,0,0,'','','',NaN,NaN,NaN,'-','none','b',2);
0690    [ax] = graph1D_jd(theta_fit,dxdrho_fit1,0,0,'','','',NaN,NaN,NaN,'-','none','c',1);
0691    [ax] = graph1D_jd(equil_in.theta(1:2:end),dxdrho(1:nr:end,1:2:end),0,0,'\theta','dx/d\rho',[strrep(equil_in.id,'_','-'),', blue: ppval, cyan: ppval-yp, red: numerical'],NaN,NaN,NaN,'--','non','r',1);   
0692    %
0693    figure('Name','d2xdthetadrho'),
0694    [ax] = graph1D_jd(theta_fit,d2xdthetadrho_fit,0,0,'','','',NaN,NaN,NaN,'-','none','b',2);
0695    [ax] = graph1D_jd(theta_fit,d2xdthetadrho_fit1,0,0,'','','',NaN,NaN,NaN,'-','none','c',1);
0696    [ax] = graph1D_jd(equil_in.theta(1:2:end),d2xdthetadrho(1:nr:end,1:2:end),0,0,'\theta','d2xdthetadrho',[strrep(equil_in.id,'_','-'),', blue: ppval, cyan: ppval-yp, red: numerical'],NaN,NaN,NaN,'--','non','r',1);   
0697    %
0698    figure('Name','d2x/drho2'),
0699    [ax] = graph1D_jd(theta_fit,d2xdrho2_fit,0,0,'','','',NaN,NaN,NaN,'-','none','b',2);
0700    [ax] = graph1D_jd(theta_fit,d2xdrho2_fit1,0,0,'','','',NaN,NaN,NaN,'-','none','c',1);
0701    [ax] = graph1D_jd(equil_in.theta(1:2:end),d2xdrho2(1:nr:end,1:2:end),0,0,'\theta','d2x/d\rho2',[strrep(equil_in.id,'_','-'),', blue: ppval, cyan: ppval-yp, red: numerical'],NaN,NaN,NaN,'--','non','r',1);   
0702    %
0703    figure('Name','y (vertical) coordinate'),
0704    [ax] = graph1D_jd(theta_fit,y_fit,0,0,'','','',NaN,NaN,NaN,'-','none','b',2);
0705    [ax] = graph1D_jd(theta_fit,y_fit1,0,0,'','','',NaN,NaN,NaN,'-','none','c',1);
0706    [ax] = graph1D_jd(equil_in.theta(1:2:end),y(1:nr:end,1:2:end),0,0,'\theta (rd)','y (m)',[strrep(equil_in.id,'_','-'),', blue: ppval, cyan: ppval-yp, red: numerical'],NaN,NaN,NaN,'--','non','r',1);
0707    %
0708    figure('Name','dy/dtheta'),
0709    [ax] = graph1D_jd(theta_fit,dydtheta_fit,0,0,'','','',NaN,NaN,NaN,'-','none','b',2);
0710    [ax] = graph1D_jd(theta_fit,dydtheta_fit1,0,0,'','','',NaN,NaN,NaN,'-','none','c',1);
0711    [ax] = graph1D_jd(equil_in.theta(1:2:end),dydtheta(1:nr:end,1:2:end),0,0,'\theta (rd)','dy/d\theta',[strrep(equil_in.id,'_','-'),', blue: ppval, cyan: ppval-yp, red: numerical'],NaN,NaN,NaN,'--','non','r',1);
0712    %
0713    figure('Name','dy/drho'),
0714    [ax] = graph1D_jd(theta_fit,dydrho_fit,0,0,'','','',NaN,NaN,NaN,'-','none','b',2);
0715    [ax] = graph1D_jd(theta_fit,dydrho_fit1,0,0,'','','',NaN,NaN,NaN,'-','none','c',1);
0716    [ax] = graph1D_jd(equil_in.theta(1:2:end),dydrho(1:nr:end,1:2:end),0,0,'\rho','dy/d\rho',[strrep(equil_in.id,'_','-'),', blue: ppval, cyan: ppval-yp, red: numerical'],NaN,NaN,NaN,'--','non','r',1);
0717    %
0718    figure('Name','d2ydthetadrho'),
0719    [ax] = graph1D_jd(theta_fit,d2ydthetadrho_fit,0,0,'','','',NaN,NaN,NaN,'-','none','b',2);
0720    [ax] = graph1D_jd(theta_fit,d2ydthetadrho_fit1,0,0,'','','',NaN,NaN,NaN,'-','none','c',1);
0721    [ax] = graph1D_jd(equil_in.theta(1:2:end),d2ydthetadrho(1:nr:end,1:2:end),0,0,'\theta','d2ydthetadrho',[strrep(equil_in.id,'_','-'),', blue: ppval, cyan: ppval-yp, red: numerical'],NaN,NaN,NaN,'--','non','r',1);   
0722    %
0723    figure('Name','d2y/drho2'),
0724    [ax] = graph1D_jd(theta_fit,d2ydrho2_fit,0,0,'','','',NaN,NaN,NaN,'-','none','b',2);
0725    [ax] = graph1D_jd(theta_fit,d2ydrho2_fit1,0,0,'','','',NaN,NaN,NaN,'-','none','c',1);
0726    [ax] = graph1D_jd(equil_in.theta(1:2:end),d2ydrho2(1:nr:end,1:2:end),0,0,'\theta','d2y/d\rho2',[strrep(equil_in.id,'_','-'),', blue: ppval, cyan: ppval-yp, red: numerical'],NaN,NaN,NaN,'--','non','r',1);   
0727    %
0728    % Magnetic equilibrium quantities
0729    %
0730    figure('Name','r position'),
0731    [ax] = graph1D_jd(theta_fit,r_fit,0,0,'','','',NaN,NaN,NaN,'-','none','b',2);
0732    [ax] = graph1D_jd(theta_fit,r_fit1,0,0,'','','',NaN,NaN,NaN,'-','none','c',1);
0733    [ax] = graph1D_jd(equil_in.theta(1:2:end),r(1:nr:end,1:2:end),0,0,'\theta (rd)','r (m)',[strrep(equil_in.id,'_','-'),', blue: ppval, cyan: ppval-yp, red: numerical'],NaN,NaN,NaN,'--','non','r',1);
0734    %
0735    figure('Name','dr/dtheta'),
0736    [ax] = graph1D_jd(theta_fit,drdtheta_fit,0,0,'','','',NaN,NaN,NaN,'-','none','b',2);
0737    [ax] = graph1D_jd(theta_fit,drdtheta_fit1,0,0,'','','',NaN,NaN,NaN,'-','none','c',1);
0738    [ax] = graph1D_jd(equil_in.theta(1:2:end),drdtheta(1:nr:end,1:2:end),0,0,'\theta (rd)','dr/d\theta',[strrep(equil_in.id,'_','-'),', blue: ppval, cyan: ppval-yp, red: numerical'],NaN,NaN,NaN,'--','non','r',1);
0739    %
0740    figure('Name','dr/drho'),
0741    [ax] = graph1D_jd(theta_fit,drdrho_fit,0,0,'','','',NaN,NaN,NaN,'-','none','b',2);
0742    [ax] = graph1D_jd(theta_fit,drdrho_fit1,0,0,'','','',NaN,NaN,NaN,'-','none','c',1);
0743    [ax] = graph1D_jd(equil_in.theta(1:2:end),drdrho(1:nr:end,1:2:end),0,0,'\rho','dr/d\rho',[strrep(equil_in.id,'_','-'),', blue: ppval, cyan: ppval-yp, red: numerical'],NaN,NaN,NaN,'--','non','r',1);
0744    %
0745    figure('Name','d2rdthetadrho'),
0746    [ax] = graph1D_jd(theta_fit,d2ydthetadrho_fit,0,0,'','','',NaN,NaN,NaN,'-','none','b',2);
0747    [ax] = graph1D_jd(theta_fit,d2ydthetadrho_fit1,0,0,'','','',NaN,NaN,NaN,'-','none','c',1);
0748    [ax] = graph1D_jd(equil_in.theta(1:2:end),d2ydthetadrho(1:nr:end,1:2:end),0,0,'\theta','d2rdthetadrho',[strrep(equil_in.id,'_','-'),', blue: ppval, cyan: ppval-yp, red: numerical'],NaN,NaN,NaN,'--','non','r',1);   
0749    %
0750    figure('Name','d2r/drho2'),
0751    [ax] = graph1D_jd(theta_fit,d2rdrho2_fit,0,0,'','','',NaN,NaN,NaN,'-','none','b',2);
0752    [ax] = graph1D_jd(theta_fit,d2rdrho2_fit1,0,0,'','','',NaN,NaN,NaN,'-','none','c',1);
0753    [ax] = graph1D_jd(equil_in.theta(1:2:end),d2rdrho2(1:nr:end,1:2:end),0,0,'\theta','d2r/d\rho2',[strrep(equil_in.id,'_','-'),', blue: ppval, cyan: ppval-yp, red: numerical'],NaN,NaN,NaN,'--','non','r',1);   
0754    %
0755    figure('Name','cos(alpha)'),
0756    [ax] = graph1D_jd(theta_fit,calpha_fit,0,0,'','','',NaN,NaN,NaN,'-','none','b',2);
0757    [ax] = graph1D_jd(theta_fit,calpha_fit1,0,0,'','','',NaN,NaN,NaN,'-','none','c',1);
0758    [ax] = graph1D_jd(equil_in.theta(1:2:end),calpha(1:nr:end,1:2:end),0,0,'\theta (rd)','cos(\alpha)',[strrep(equil_in.id,'_','-'),', blue: ppval, cyan: ppval-yp, red: numerical'],NaN,NaN,NaN,'--','non','r',1);
0759    %
0760    figure('Name','dcos(alpha)/dtheta'),
0761    [ax] = graph1D_jd(theta_fit,dcalphadtheta_fit,0,0,'','','',NaN,NaN,NaN,'-','none','b',2);
0762    [ax] = graph1D_jd(theta_fit,dcalphadtheta_fit1,0,0,'','','',NaN,NaN,NaN,'-','none','c',1);
0763    [ax] = graph1D_jd(equil_in.theta(1:2:end),dcalphadtheta(1:nr:end,1:2:end),0,0,'\theta (rd)','dcos(\alpha)/d\theta',[strrep(equil_in.id,'_','-'),', blue: ppval, cyan: ppval-yp, red: numerical'],NaN,NaN,NaN,'--','non','r',1);
0764    %
0765    figure('Name','dcos(alpha)/drho'),
0766    [ax] = graph1D_jd(theta_fit,dcalphadrho_fit,0,0,'','','',NaN,NaN,NaN,'-','none','b',2);
0767    [ax] = graph1D_jd(theta_fit,dcalphadrho_fit1,0,0,'','','',NaN,NaN,NaN,'-','none','c',1);
0768    [ax] = graph1D_jd(equil_in.theta(1:2:end),dcalphadrho(1:nr:end,1:2:end),0,0,'\rho','dcos(\alpha)/d\rho',[strrep(equil_in.id,'_','-'),', blue: ppval, cyan: ppval-yp, red: numerical'],NaN,NaN,NaN,'--','none','r',1);
0769    %
0770    figure('Name','d2cos(alpha)/dtheta/drho'),
0771    [ax] = graph1D_jd(theta_fit,d2calphadthetadrho_fit,0,0,'','','',NaN,NaN,NaN,'-','none','b',2);
0772    [ax] = graph1D_jd(theta_fit,d2calphadthetadrho_fit1,0,0,'','','',NaN,NaN,NaN,'-','none','c',1);
0773    [ax] = graph1D_jd(equil_in.theta(1:2:end),d2calphadthetadrho(1:nr:end,1:2:end),0,0,'\theta','d2cos(alpha)/d\theta/d\rho',[strrep(equil_in.id,'_','-'),', blue: ppval, cyan: ppval-yp, red: numerical'],NaN,NaN,NaN,'--','non','r',1);   
0774    %
0775    figure('Name','d2cos(alpha)/drho2'),
0776    [ax] = graph1D_jd(theta_fit,d2calphadrho2_fit,0,0,'','','',NaN,NaN,NaN,'-','none','b',2);
0777    [ax] = graph1D_jd(theta_fit,d2calphadrho2_fit1,0,0,'','','',NaN,NaN,NaN,'-','none','c',1);
0778    [ax] = graph1D_jd(equil_in.theta(1:2:end),d2calphadrho2(1:nr:end,1:2:end),0,0,'\theta','d2cos(alpha)/d\rho2',[strrep(equil_in.id,'_','-'),', blue: ppval, cyan: ppval-yp, red: numerical'],NaN,NaN,NaN,'--','non','r',1);   
0779    %
0780    figure('Name','sin(alpha)'),
0781    [ax] = graph1D_jd(theta_fit,salpha_fit,0,0,'','','',NaN,NaN,NaN,'-','none','b',2);
0782    [ax] = graph1D_jd(theta_fit,salpha_fit1,0,0,'','','',NaN,NaN,NaN,'-','none','c',1);
0783    [ax] = graph1D_jd(equil_in.theta(1:2:end),salpha(1:nr:end,1:2:end),0,0,'\theta (rd)','sin(\alpha)',[strrep(equil_in.id,'_','-'),', blue: ppval, cyan: ppval-yp, red: numerical'],NaN,NaN,NaN,'--','non','r',1);
0784    %
0785    figure('Name','dsin(alpha)/dtheta'),
0786    [ax] = graph1D_jd(theta_fit,dsalphadtheta_fit,0,0,'','','',NaN,NaN,NaN,'-','none','b',2);
0787    [ax] = graph1D_jd(theta_fit,dsalphadtheta_fit1,0,0,'','','',NaN,NaN,NaN,'-','none','c',1);
0788    [ax] = graph1D_jd(equil_in.theta(1:2:end),dsalphadtheta(1:nr:end,1:2:end),0,0,'\theta (rd)','dsin(\alpha)/d\theta',[strrep(equil_in.id,'_','-'),', blue: ppval, cyan: ppval-yp, red: numerical'],NaN,NaN,NaN,'--','non','r',1);
0789    %
0790    figure('Name','dsin(alpha)/drho'),
0791    [ax] = graph1D_jd(theta_fit,dsalphadrho_fit,0,0,'','','',NaN,NaN,NaN,'-','none','b',2);
0792    [ax] = graph1D_jd(theta_fit,dsalphadrho_fit1,0,0,'','','',NaN,NaN,NaN,'-','none','c',1);
0793    [ax] = graph1D_jd(equil_in.theta(1:2:end),dsalphadrho(1:nr:end,1:2:end),0,0,'\rho','dsin(\alpha)/d\rho',[strrep(equil_in.id,'_','-'),', blue: ppval, cyan: ppval-yp, red: numerical'],NaN,NaN,NaN,'--','non','r',1);
0794    %
0795    figure('Name','d2sin(alpha)/dtheta/drho'),
0796    [ax] = graph1D_jd(theta_fit,d2salphadthetadrho_fit,0,0,'','','',NaN,NaN,NaN,'-','none','b',2);
0797    [ax] = graph1D_jd(theta_fit,d2salphadthetadrho_fit1,0,0,'','','',NaN,NaN,NaN,'-','none','c',1);
0798    [ax] = graph1D_jd(equil_in.theta(1:2:end),d2salphadthetadrho(1:nr:end,1:2:end),0,0,'\theta','dsin(\alpha)/d\theta/d\rho',[strrep(equil_in.id,'_','-'),', blue: ppval, cyan: ppval-yp, red: numerical'],NaN,NaN,NaN,'--','non','r',1);   
0799    %
0800    figure('Name','d2sin(alpha)/drho2'),
0801    [ax] = graph1D_jd(theta_fit,d2salphadrho2_fit,0,0,'','','',NaN,NaN,NaN,'-','none','b',2);
0802    [ax] = graph1D_jd(theta_fit,d2salphadrho2_fit1,0,0,'','','',NaN,NaN,NaN,'-','none','c',1);
0803    [ax] = graph1D_jd(equil_in.theta(1:2:end),d2salphadrho2(1:nr:end,1:2:end),0,0,'\theta','d2sin(\alpha)/d\rho2',[strrep(equil_in.id,'_','-'),', blue: ppval, cyan: ppval-yp, red: numerical'],NaN,NaN,NaN,'--','non','r',1);   
0804    %
0805    figure('Name','cos^2(alpha)+sin^2(alpha)'),
0806    [ax] = graph1D_jd(theta_fit,salpha_fit.^2 + calpha_fit.^2,0,0,'','','',NaN,NaN,NaN,'-','none','b',2);
0807    [ax] = graph1D_jd(equil_in.theta(1:2:end),salpha(1:nr:end,1:2:end).^2 + calpha(1:nr:end,1:2:end).^2,0,0,'\theta (rd)','cos^2(\alpha)+sin^2(\alpha)',[strrep(equil_in.id,'_','-'),', blue: ppval, cyan: ppval-yp, red: numerical'],NaN,NaN,NaN,'--','non','r',1);
0808    %
0809    figure('Name','\nabla*\rho'),
0810    [ax] = graph1D_jd(theta_fit,gradrho_fit(1:end,:),0,0,'','','',NaN,NaN,NaN,'-','none','b',2);
0811    [ax] = graph1D_jd(theta_fit,gradrho_fit1(1:end,:),0,0,'','','',NaN,NaN,NaN,'-','none','c',1);
0812    [ax] = graph1D_jd(equil_in.theta(1:2:end),gradrho(1:nr:end,1:2:end),0,0,'\theta (rd)','\nabla*\rho',[strrep(equil_in.id,'_','-'),', blue: ppval, cyan: ppval-yp, red: numerical'],NaN,NaN,NaN,'--','non','r',1);
0813    %
0814    figure('Name','d\nabla*\rho/dtheta'),
0815    [ax] = graph1D_jd(theta_fit,dgradrhodtheta_fit,0,0,'','','',NaN,NaN,NaN,'-','none','b',2);
0816    [ax] = graph1D_jd(theta_fit,dgradrhodtheta_fit1,0,0,'','','',NaN,NaN,NaN,'-','none','c',1);
0817    [ax] = graph1D_jd(equil_in.theta(1:2:end),dgradrhodtheta(1:nr:end,1:2:end),0,0,'\theta (rd)','d\nabla*\rho/d\theta',[strrep(equil_in.id,'_','-'),', blue: ppval, cyan: ppval-yp, red: numerical'],NaN,NaN,NaN,'--','non','r',1);
0818    %
0819    figure('Name','d\nabla*\rho/drho'),
0820    [ax] = graph1D_jd(theta_fit,dgradrhodrho_fit,0,0,'','','',NaN,NaN,NaN,'-','none','b',2);
0821    [ax] = graph1D_jd(theta_fit,dgradrhodrho_fit1,0,0,'','','',NaN,NaN,NaN,'-','none','c',1);
0822    [ax] = graph1D_jd(equil_in.theta(1:2:end),dgradrhodrho(1:nr:end,1:2:end),0,0,'\rho','d\nabla*\rho/d\rho',[strrep(equil_in.id,'_','-'),', blue: ppval, cyan: ppval-yp, red: numerical'],NaN,NaN,NaN,'--','non','r',1);
0823    %
0824    figure('Name','d2\nabla*\rho/dtheta/drho'),
0825    [ax] = graph1D_jd(theta_fit,d2gradrhodthetadrho_fit,0,0,'','','',NaN,NaN,NaN,'-','none','b',2);
0826    [ax] = graph1D_jd(theta_fit,d2gradrhodthetadrho_fit1,0,0,'','','',NaN,NaN,NaN,'-','none','c',1);
0827    [ax] = graph1D_jd(equil_in.theta(1:2:end),d2gradrhodthetadrho(1:nr:end,1:2:end),0,0,'\theta','d2\nabla*\rhodthetadrho',[strrep(equil_in.id,'_','-'),', blue: ppval, cyan: ppval-yp, red: numerical'],NaN,NaN,NaN,'--','non','r',1);   
0828    %
0829    figure('Name','d2\nabla*\rho/drho2'),
0830    [ax] = graph1D_jd(theta_fit,d2ydrho2_fit,0,0,'','','',NaN,NaN,NaN,'-','none','b',2);
0831    [ax] = graph1D_jd(theta_fit,d2ydrho2_fit1,0,0,'','','',NaN,NaN,NaN,'-','none','c',1);
0832    [ax] = graph1D_jd(equil_in.theta(1:2:end),d2ydrho2(1:nr:end,1:2:end),0,0,'\theta','d2\nabla*\rho/d\rho2',[strrep(equil_in.id,'_','-'),', blue: ppval, cyan: ppval-yp, red: numerical'],NaN,NaN,NaN,'--','non','r',1);   
0833    %
0834    % Magnetic fields
0835    %
0836    figure('Name','Toroidal (axial) magnetic field Bz'),
0837    [ax] = graph1D_jd(theta_fit,Bz_fit,0,0,'\rho','\rho','',NaN,NaN,NaN,'-','none','b',2);
0838    [ax] = graph1D_jd(theta_fit,Bz_fit1,0,0,'\rho','\rho','',NaN,NaN,NaN,'-','none','c',1);
0839    [ax] = graph1D_jd(equil_in.theta(1:2:end),Bz(1:nr:end,1:2:end),0,0,'\theta (rd)','B_{\phi} (T)',[strrep(equil_in.id,'_','-'),', blue: ppval, cyan: ppval-yp, red: numerical'],NaN,NaN,NaN,'--','non','r',1);
0840    %
0841    figure('Name','dBz/dtheta'),
0842    [ax] = graph1D_jd(theta_fit,dBzdtheta_fit,0,0,'','','',NaN,NaN,NaN,'-','none','b',2);
0843    [ax] = graph1D_jd(theta_fit,dBzdtheta_fit1,0,0,'','','',NaN,NaN,NaN,'-','none','c',1);
0844    [ax] = graph1D_jd(equil_in.theta(1:2:end),dBzdtheta(1:nr:end,1:2:end),0,0,'\theta (rd)','dBz/d\theta',[strrep(equil_in.id,'_','-'),', blue: ppval, cyan: ppval-yp, red: numerical'],NaN,NaN,NaN,'--','non','r',1);
0845    %
0846    figure('Name','dBz/drho'),
0847    [ax] = graph1D_jd(theta_fit,dBzdrho_fit,0,0,'','','',NaN,NaN,NaN,'-','none','b',2);
0848    [ax] = graph1D_jd(theta_fit,dBzdrho_fit1,0,0,'','','',NaN,NaN,NaN,'-','none','c',1);
0849    [ax] = graph1D_jd(equil_in.theta(1:2:end),dBzdrho(1:nr:end,1:2:end),0,0,'\psi_{n}','dBz/d\rho',[strrep(equil_in.id,'_','-'),', blue: ppval, cyan: ppval-yp, red: numerical'],NaN,NaN,NaN,'--','non','r',1);
0850    %
0851    figure('Name','d2Bz/dtheta/drho'),
0852    [ax] = graph1D_jd(theta_fit,d2Bzdthetadrho_fit,0,0,'','','',NaN,NaN,NaN,'-','none','b',2);
0853    [ax] = graph1D_jd(theta_fit,d2Bzdthetadrho_fit1,0,0,'','','',NaN,NaN,NaN,'-','none','c',1);
0854    [ax] = graph1D_jd(equil_in.theta(1:2:end),d2Bzdthetadrho(1:nr:end,1:2:end),0,0,'\theta','d2Bzdthetadrho',[strrep(equil_in.id,'_','-'),', blue: ppval, cyan: ppval-yp, red: numerical'],NaN,NaN,NaN,'--','non','r',1);   
0855    %
0856    figure('Name','d2Bz/drho2'),
0857    [ax] = graph1D_jd(theta_fit,d2Bzdrho2_fit,0,0,'','','',NaN,NaN,NaN,'-','none','b',2);
0858    [ax] = graph1D_jd(theta_fit,d2Bzdrho2_fit1,0,0,'','','',NaN,NaN,NaN,'-','none','c',1);
0859    [ax] = graph1D_jd(equil_in.theta(1:2:end),d2Bzdrho2(1:nr:end,1:2:end),0,0,'\theta','d2Bz/d\rho2',[strrep(equil_in.id,'_','-'),', blue: ppval, cyan: ppval-yp, red: numerical'],NaN,NaN,NaN,'--','non','r',1);   
0860    %
0861    figure('Name','Poloidal magnetic field BP'),
0862    [ax] = graph1D_jd(theta_fit,BP_fit,0,0,'\rho','\rho','',NaN,NaN,NaN,'-','none','b',2);
0863    [ax] = graph1D_jd(theta_fit,BP_fit1,0,0,'\rho','\rho','',NaN,NaN,NaN,'-','none','c',1);
0864    [ax] = graph1D_jd(equil_in.theta(1:2:end),BP(1:nr:end,1:2:end),0,0,'\theta (rd)','B_P (T)',[strrep(equil_in.id,'_','-'),', blue: ppval, cyan: ppval-yp, red: numerical'],NaN,NaN,NaN,'--','non','r',1);
0865    %
0866    figure('Name','dBP/dtheta'),
0867    [ax] = graph1D_jd(theta_fit,dBPdtheta_fit,0,0,'\rho','\rho','',NaN,NaN,NaN,'-','none','b',2);
0868    [ax] = graph1D_jd(theta_fit,dBPdtheta_fit1,0,0,'\rho','\rho','',NaN,NaN,NaN,'-','none','c',1);
0869    [ax] = graph1D_jd(equil_in.theta(1:2:end),dBPdtheta(1:nr:end,1:2:end),0,0,'\theta (rd)','dB_P/d\theta (T)',[strrep(equil_in.id,'_','-'),', blue: ppval, cyan: ppval-yp, red: numerical'],NaN,NaN,NaN,'--','non','r',1);
0870    %
0871    figure('Name','dBP/drho'),
0872    [ax] = graph1D_jd(theta_fit,dBPdrho_fit,0,0,'\rho','\rho','',NaN,NaN,NaN,'-','none','b',2);
0873    [ax] = graph1D_jd(theta_fit,dBPdrho_fit1,0,0,'\rho','\rho','',NaN,NaN,NaN,'-','none','c',1);
0874    [ax] = graph1D_jd(equil_in.theta(1:2:end),dBPdrho(1:nr:end,1:2:end),0,0,'\theta (rd)','dB_P/d\rho (T)',[strrep(equil_in.id,'_','-'),', blue: ppval, cyan: ppval-yp, red: numerical'],NaN,NaN,NaN,'--','non','r',1);
0875    %
0876    figure('Name','d2BP/dtheta/drho'),
0877    [ax] = graph1D_jd(theta_fit,d2BPdthetadrho_fit,0,0,'','','',NaN,NaN,NaN,'-','none','b',2);
0878    [ax] = graph1D_jd(theta_fit,d2BPdthetadrho_fit1,0,0,'','','',NaN,NaN,NaN,'-','none','c',1);
0879    [ax] = graph1D_jd(equil_in.theta(1:2:end),d2BPdthetadrho(1:nr:end,1:2:end),0,0,'\theta','d2BPdthetadrho',[strrep(equil_in.id,'_','-'),', blue: ppval, cyan: ppval-yp, red: numerical'],NaN,NaN,NaN,'--','non','r',1);   
0880    %
0881    figure('Name','d2BP/drho2'),
0882    [ax] = graph1D_jd(theta_fit,d2BPdrho2_fit,0,0,'','','',NaN,NaN,NaN,'-','none','b',2);
0883    [ax] = graph1D_jd(theta_fit,d2BPdrho2_fit1,0,0,'','','',NaN,NaN,NaN,'-','none','c',1);
0884    [ax] = graph1D_jd(equil_in.theta(1:2:end),d2BPdrho2(1:nr:end,1:2:end),0,0,'\theta','d2BP/d\rho2',[strrep(equil_in.id,'_','-'),', blue: ppval, cyan: ppval-yp, red: numerical'],NaN,NaN,NaN,'--','non','r',1);   
0885    %
0886    figure('Name','Total magnetic field'),
0887    [ax] = graph1D_jd(theta_fit,B_fit,0,0,'\rho','\rho','',NaN,NaN,NaN,'-','none','b',2);
0888    [ax] = graph1D_jd(theta_fit,B_fit1,0,0,'\rho','\rho','',NaN,NaN,NaN,'-','none','c',1);
0889    [ax] = graph1D_jd(equil_in.theta(1:2:end),B(1:nr:end,1:2:end),0,0,'\theta (rd)','B (T)',[strrep(equil_in.id,'_','-'),', blue: ppval, cyan: ppval-yp, red: numerical'],NaN,NaN,NaN,'--','non','r',1);
0890    %
0891    figure('Name','dB/dtheta'),
0892    [ax] = graph1D_jd(theta_fit,dBdtheta_fit,0,0,'\rho','\rho','',NaN,NaN,NaN,'-','none','b',2);
0893    [ax] = graph1D_jd(theta_fit,dBdtheta_fit1,0,0,'\rho','\rho','',NaN,NaN,NaN,'-','none','c',1);
0894    [ax] = graph1D_jd(equil_in.theta(1:2:end),dBdtheta(1:nr:end,1:2:end),0,0,'\theta (rd)','dB/d\theta (T)',[strrep(equil_in.id,'_','-'),', blue: ppval, cyan: ppval-yp, red: numerical'],NaN,NaN,NaN,'--','non','r',1);
0895    %
0896    figure('Name','dB/drho'),
0897    [ax] = graph1D_jd(theta_fit,dBdrho_fit,0,0,'\rho','\rho','',NaN,NaN,NaN,'-','none','b',2);
0898    [ax] = graph1D_jd(theta_fit,dBdrho_fit1,0,0,'\rho','\rho','',NaN,NaN,NaN,'-','none','c',1);
0899    [ax] = graph1D_jd(equil_in.theta(1:2:end),dBdrho(1:nr:end,1:2:end),0,0,'\theta (rd)','dB/d\rho (T)',[strrep(equil_in.id,'_','-'),', blue: ppval, cyan: ppval-yp, red: numerical'],NaN,NaN,NaN,'--','non','r',1);
0900    %
0901    figure('Name','d2B/dtheta/drho'),
0902    [ax] = graph1D_jd(theta_fit,d2Bdthetadrho_fit,0,0,'','','',NaN,NaN,NaN,'-','none','b',2);
0903    [ax] = graph1D_jd(theta_fit,d2Bdthetadrho_fit1,0,0,'','','',NaN,NaN,NaN,'-','none','c',1);
0904    [ax] = graph1D_jd(equil_in.theta(1:2:end),d2Bdthetadrho(1:nr:end,1:2:end),0,0,'\theta','d2Bdthetadrho',[strrep(equil_in.id,'_','-'),', blue: ppval, cyan: ppval-yp, red: numerical'],NaN,NaN,NaN,'--','non','r',1);   
0905    %
0906    figure('Name','d2B/drho2'),
0907    [ax] = graph1D_jd(theta_fit,d2Bdrho2_fit,0,0,'','','',NaN,NaN,NaN,'-','none','b',2);
0908    [ax] = graph1D_jd(theta_fit,d2Bdrho2_fit1,0,0,'','','',NaN,NaN,NaN,'-','none','c',1);
0909    [ax] = graph1D_jd(equil_in.theta(1:2:end),d2Bdrho2(1:nr:end,1:2:end),0,0,'\theta','d2B/d\rho2',[strrep(equil_in.id,'_','-'),', blue: ppval, cyan: ppval-yp, red: numerical'],NaN,NaN,NaN,'--','non','r',1);   
0910    %
0911    xmin = min(min(x));
0912    xmax = max(max(x));
0913    ymin = min(min(y));
0914    ymax = max(max(y));
0915    %
0916    figure('Name','Magnetic equilibrium flux surface contour')
0917    %
0918    [ax] = graph1D_jd(x(1:nr:end,1:2:end),y(1:nr:end,1:2:end),0,0,'','','',NaN,[xmin,xmax],[ymin,ymax],'-','none','r',2);drawnow
0919    [ax] = graph1D_jd(x(1:nr:end,:)',y(1:nr:end,:)',0,0,'x(m)','y(m)','Flux-surfaces labeled by \rho=0,0.1,0.2...1',NaN,[xmin,xmax],[ymin,ymax],'-','none','r',2);
0920    %
0921    [ax] = graph1D_jd(x_fit(1:nr:end,1:2:end),y_fit(1:nr:end,1:2:end),0,0,'','','',NaN,[xmin,xmax],[ymin,ymax],'-','none','b',1);drawnow
0922    [ax] = graph1D_jd(x_fit(:,:)',y_fit(:,:)',0,0,'x(m)','y(m)','Flux-surfaces labeled by \rho=0,0.1,0.2...1',NaN,[xmin,xmax],[ymin,ymax],'-','none','b',1);
0923    %
0924    [ax] = graph1D_jd(x_fit1(1:nr:end,1:2:end),y_fit1(1:nr:end,1:2:end),0,0,'','','',NaN,[xmin,xmax],[ymin,ymax],'--','none','c',1);drawnow
0925    [ax] = graph1D_jd(x_fit1(:,:)',y_fit1(:,:)',0,0,'x(m)','y(m)','Flux-surfaces labeled by \rho=0,0.1,0.2...1',NaN,[xmin,xmax],[ymin,ymax],'--','none','c',1);
0926    axis('equal');drawnow  
0927    %
0928 else
0929     equilval = equilval_yp(equil_fit,rho,theta);
0930     %
0931     Te = equilval.Te;
0932     dTedrho = equilval.dTedrho;
0933     d2Tedrho2 = equilval.d2Tedrho2;
0934     %
0935     ne = equilval.ne;
0936     dnedrho = equilval.dnedrho;
0937     d2nedrho2 = equilval.d2nedrho2;
0938     %
0939     zTi = equilval.zTi;
0940     dzTidrho = equilval.dzTidrho;
0941     d2zTidrho2 = equilval.d2zTidrho2;
0942     %
0943     zni = equilval.zni;
0944     dznidrho = equilval.dznidrho;
0945     d2znidrho2 = equilval.d2znidrho2;
0946     %
0947     x = equilval.x;
0948     dxdtheta = equilval.dxdtheta;
0949     d2xdtheta2 = equilval.d2xdtheta2;
0950     dxdrho = equilval.dxdrho;
0951     d2xdthetadrho = equilval.d2xdthetadrho;
0952     d2xdrho2 = equilval.d2xdrho2;
0953     %
0954     y = equilval.y;
0955     dydtheta = equilval.dydtheta;
0956     d2ydtheta2 = equilval.d2ydtheta2;
0957     dydrho = equilval.dydrho;
0958     d2ydthetadrho = equilval.d2ydthetadrho;
0959     d2ydrho2 = equilval.d2ydrho2;
0960     %
0961     r = equilval.r;
0962     drdtheta = equilval.drdtheta;
0963     d2rdtheta2 = equilval.d2rdtheta2;
0964     drdrho = equilval.drdrho;
0965     d2rdthetadrho = equilval.d2rdthetadrho;
0966     d2rdrho2 = equilval.d2rdrho2;
0967     %
0968     salpha = equilval.salpha;
0969     dsalphadtheta = equilval.dsalphadtheta;
0970     d2salphadtheta2 = equilval.d2salphadtheta2;
0971     dsalphadrho = equilval.dsalphadrho;
0972     d2salphadthetadrho = equilval.d2salphadthetadrho;
0973     d2salphadrho2 = equilval.d2salphadrho2;
0974     %
0975     calpha = equilval.calpha;
0976     dcalphadtheta = equilval.dcalphadtheta;
0977     d2calphadtheta2 = equilval.d2calphadtheta2;
0978     dcalphadrho = equilval.dcalphadrho;
0979     d2calphadthetadrho = equilval.d2calphadthetadrho;
0980     d2calphadrho2 = equilval.d2calphadrho2;
0981     %
0982     Bx = equilval.Bx;
0983     dBxdtheta = equilval.dBxdtheta;
0984     d2Bxdtheta2 = equilval.d2Bxdtheta2;
0985     dBxdrho = equilval.dBxdrho;
0986     d2Bxdthetadrho = equilval.d2Bxdthetadrho;
0987     d2Bxdrho2 = equilval.d2Bxdrho2;
0988     %
0989     By = equilval.By;
0990     dBydtheta = equilval.dBydtheta;
0991     d2Bydtheta2 = equilval.d2Bydtheta2;
0992     dBydrho = equilval.dBydrho;
0993     d2Bydthetadrho = equilval.d2Bydthetadrho;
0994     d2Bydrho2 = equilval.d2Bydrho2;
0995     %
0996     Bz = equilval.Bz;
0997     dBzdtheta = equilval.dBzdtheta;
0998     d2Bzdtheta2 = equilval.d2Bzdtheta2;
0999     dBzdrho = equilval.dBzdrho;
1000     d2Bzdthetadrho = equilval.d2Bzdthetadrho;
1001     d2Bzdrho2 = equilval.d2Bzdrho2;
1002     %
1003     BP = equilval.BP;
1004     dBPdtheta = equilval.dBPdtheta;
1005     d2BPdtheta2 = equilval.d2BPdtheta2;
1006     dBPdrho = equilval.dBPdrho;
1007     d2BPdthetadrho = equilval.d2BPdthetadrho;
1008     d2BPdrho2 = equilval.d2BPdrho2;
1009     %
1010     B = equilval.B;
1011     dBdtheta = equilval.dBdtheta;
1012     d2Bdtheta2 = equilval.d2Bdtheta2;
1013     dBdrho = equilval.dBdrho;
1014     d2Bdthetadrho = equilval.d2Bdthetadrho;
1015     d2Bdrho2 = equilval.d2Bdrho2;
1016     %
1017     gradrho = equilval.gradrho;
1018     dgradrhodtheta = equilval.dgradrhodtheta;
1019     d2gradrhodtheta2 = equilval.d2gradrhodtheta2;
1020     dgradrhodrho = equilval.dgradrhodrho;
1021     d2gradrhodthetadrho = equilval.d2gradrhodthetadrho;
1022     d2gradrhodrho2 = equilval.d2gradrhodrho2;
1023     %
1024     Upsilon = equilval.Upsilon;
1025     dUpsilondtheta = equilval.dUpsilondtheta;
1026     d2Upsilondtheta2 = equilval.d2Upsilondtheta2;
1027     dUpsilondrho = equilval.dUpsilondrho;
1028     d2Upsilondthetadrho = equilval.d2Upsilondthetadrho;
1029     d2Upsilondrho2 = equilval.d2Upsilondrho2;
1030     %
1031     %psin = equilval.psin;
1032     %dpsindrho = equilval.dpsindrho;
1033     %dxstardpsin = equilval.dxstardpsin;
1034     %ddxstardpsindrho = equilval.ddxstardpsindrho;
1035     %
1036     if nargout ~= 2,
1037         disp('----------------------------');
1038         if isfield(equil_fit,'Te_fit'),
1039             disp(['equil_fit.Te_fit.pp_f.coefs(1) = ',num2str(equil_fit.Te_fit.pp_f.coefs(1))]);
1040             disp(['equil_fit.Te_fit.pp_dfdrho.coefs(1) = ',num2str(equil_fit.Te_fit.pp_dfdrho.coefs(1))]);
1041             disp(['equil_fit.Te_fit.pp_d2fdrho2.coefs(1) = ',num2str(equil_fit.Te_fit.pp_d2fdrho2.coefs(1))]);
1042         end
1043         disp(['equil_fit.Bz_fit.n = ',num2str(equil_fit.Bz_fit.n(end))]);
1044         disp(['equil_fit.Bz_fit.bn.coefs(1) = ',num2str(equil_fit.Bz_fit.pp_bn.coefs(1))]);
1045         disp(['equil_fit.Bz_fit.dbndrho.coefs(1) = ',num2str(equil_fit.Bz_fit.pp_dbndrho.coefs(1))]);
1046         disp(['equil_fit.Bz_fit.d2bndrho2.coefs(1) = ',num2str(equil_fit.Bz_fit.pp_d2bndrho2.coefs(1))]);
1047         disp('----------------------------');
1048         if isfield(equil_fit,'Te_fit'),
1049             disp(['Te(',num2str(rho),'),dTedrho(',num2str(rho),'),d2Tedrho2(',num2str(rho),') (keV) = ',num2str(Te),',',num2str(dTedrho),',',num2str(d2Tedrho2)]);
1050         end
1051         if isfield(equil_fit,'ne_fit'),
1052             disp(['ne(',num2str(rho),'),dnedrho(',num2str(rho),'),d2nedrho2(',num2str(rho),') (10+19 m-3) = ',num2str(ne/1e19),',',num2str(dnedrho/1e19),',',num2str(d2nedrho2/1e19)]);
1053         end
1054         if isfield(equil_fit,'zTi_fit'),
1055             for is = 1:ns,
1056                 disp(['zTi(',int2str(is),',',num2str(rho),'),dzTidrho(',int2str(is),',',num2str(rho),'),d2zTidrho2(',int2str(is),' ,',num2str(rho),') (keV) = ',num2str(zTi(is)),',',num2str(dzTidrho(is)),',',num2str(d2zTidrho2(is))]);
1057             end
1058         end
1059         if isfield(equil_fit,'zni_fit'),
1060             for is = 1:ns,
1061                 disp(['zni(',int2str(is),',',num2str(rho),'),dznidrho(',int2str(is),',',num2str(rho),'),d2znidrho2(',int2str(is),' ,',num2str(rho),') (10+19 m-3) = ',num2str(zni(is)/1e19),',',num2str(dznidrho(is)/1e19),',',num2str(d2znidrho2(is)/1e19)]);
1062             end
1063         end
1064         %
1065         disp('----------------------------');
1066         disp(['x(',num2str(rho),',',num2str(theta),'),dxdrho(',num2str(rho),',',num2str(theta),'),dxdtheta(',num2str(rho),',',num2str(theta),') (m) = ',num2str(x),',',num2str(dxdrho),',',num2str(dxdtheta)]);
1067         disp(['d2xdtheta2(',num2str(rho),',',num2str(theta),'),d2xdthetadrho(',num2str(rho),',',num2str(theta),'),d2xdrho2(',num2str(rho),',',num2str(theta),') (m) = ',num2str(d2xdtheta2),',',num2str(d2xdthetadrho),',',num2str(d2xdrho2)]);
1068         disp('----------------------------');
1069         disp(['y(',num2str(rho),',',num2str(theta),'),dydrho(',num2str(rho),',',num2str(theta),'),dydtheta(',num2str(rho),',',num2str(theta),') (m) = ',num2str(y),',',num2str(dydrho),',',num2str(dydtheta)]);
1070         disp(['d2ydtheta2(',num2str(rho),',',num2str(theta),'),d2ydthetadrho(',num2str(rho),',',num2str(theta),'),d2ydrho2(',num2str(rho),',',num2str(theta),') (m) = ',num2str(d2ydtheta2),',',num2str(d2ydthetadrho),',',num2str(d2ydrho2)]);
1071         disp('----------------------------');
1072         disp(['r(',num2str(rho),',',num2str(theta),'),drdrho(',num2str(rho),',',num2str(theta),'),drdtheta(',num2str(rho),',',num2str(theta),') (m) = ',num2str(r),',',num2str(drdrho),',',num2str(drdtheta)]);
1073         disp(['d2rdtheta2(',num2str(rho),',',num2str(theta),'),d2rdthetadrho(',num2str(rho),',',num2str(theta),'),d2rdrho2(',num2str(rho),',',num2str(theta),') (m) = ',num2str(d2rdtheta2),',',num2str(d2rdthetadrho),',',num2str(d2rdrho2)]);
1074         disp('----------------------------');
1075         disp(['Bx(',num2str(rho),',',num2str(theta),'),dBxdrho(',num2str(rho),',',num2str(theta),'),dBxdtheta(',num2str(rho),',',num2str(theta),') (T) = ',num2str(Bx),',',num2str(dBxdrho),',',num2str(dBxdtheta)]);
1076         disp(['d2Bxdtheta2(',num2str(rho),',',num2str(theta),'),d2Bxdthetadrho(',num2str(rho),',',num2str(theta),'),d2Bxdrho2(',num2str(rho),',',num2str(theta),')  (T) = ',num2str(d2Bxdtheta2),',',num2str(d2Bxdthetadrho),',',num2str(d2Bxdrho2)]);
1077         disp('----------------------------');
1078         disp(['By(',num2str(rho),',',num2str(theta),'),dBydrho(',num2str(rho),',',num2str(theta),'),dBydtheta(',num2str(rho),',',num2str(theta),') (T) = ',num2str(By),',',num2str(dBydrho),',',num2str(dBydtheta)]);
1079         disp(['d2Bydtheta2(',num2str(rho),',',num2str(theta),'),d2Bydthetadrho(',num2str(rho),',',num2str(theta),'),d2Bydrho2(',num2str(rho),',',num2str(theta),')  (T) = ',num2str(d2Bydtheta2),',',num2str(d2Bydthetadrho),',',num2str(d2Bydrho2)]);
1080         disp('----------------------------');
1081         disp(['Bz(',num2str(rho),',',num2str(theta),'),dBzdrho(',num2str(rho),',',num2str(theta),'),dBzdtheta(',num2str(rho),',',num2str(theta),') (T) = ',num2str(Bz),',',num2str(dBzdrho),',',num2str(dBzdtheta)]);
1082         disp(['d2Bzdtheta2(',num2str(rho),',',num2str(theta),'),d2Bzdthetadrho(',num2str(rho),',',num2str(theta),'),d2Bzdrho2(',num2str(rho),',',num2str(theta),')  (T) = ',num2str(d2Bzdtheta2),',',num2str(d2Bzdthetadrho),',',num2str(d2Bzdrho2)]);
1083         disp('----------------------------');
1084         disp(['BP(',num2str(rho),',',num2str(theta),'),dBPdrho(',num2str(rho),',',num2str(theta),'),dBPdtheta(',num2str(rho),',',num2str(theta),') (T) = ',num2str(BP),',',num2str(dBPdrho),',',num2str(dBPdtheta)]);
1085         disp(['d2BPdtheta2(',num2str(rho),',',num2str(theta),'),d2BPdthetadrho(',num2str(rho),',',num2str(theta),'),d2BPdrho2(',num2str(rho),',',num2str(theta),')  (T) = ',num2str(d2BPdtheta2),',',num2str(d2BPdthetadrho),',',num2str(d2BPdrho2)]);
1086         disp('----------------------------');
1087         disp(['B(',num2str(rho),',',num2str(theta),'),dBdrho(',num2str(rho),',',num2str(theta),'),dBdtheta(',num2str(rho),',',num2str(theta),') (T) = ',num2str(B),',',num2str(dBdrho),',',num2str(dBdtheta)]);
1088         disp(['d2Bdtheta2(',num2str(rho),',',num2str(theta),'),d2Bdthetadrho(',num2str(rho),',',num2str(theta),'),d2Bdrho2(',num2str(rho),',',num2str(theta),')  (T) = ',num2str(d2Bdtheta2),',',num2str(d2Bdthetadrho),',',num2str(d2Bdrho2)]);
1089         disp('----------------------------');
1090         disp(['calpha(',num2str(rho),',',num2str(theta),'),dcalphadrho(',num2str(rho),',',num2str(theta),'),dcalphadtheta(',num2str(rho),',',num2str(theta),') = ',num2str(calpha),',',num2str(dcalphadrho),',',num2str(dcalphadtheta)]);
1091         disp(['d2calphadtheta2(',num2str(rho),',',num2str(theta),'),d2calphadthetadrho(',num2str(rho),',',num2str(theta),'),d2calphadrho2(',num2str(rho),',',num2str(theta),') = ',num2str(d2calphadtheta2),',',num2str(d2calphadthetadrho),',',num2str(d2calphadrho2)]);
1092         disp(['salpha(',num2str(rho),',',num2str(theta),'),dsalphadrho(',num2str(rho),',',num2str(theta),'),dsalphadtheta(',num2str(rho),',',num2str(theta),') = ',num2str(salpha),',',num2str(dsalphadrho),',',num2str(dsalphadtheta)]);
1093         disp(['d2salphadtheta2(',num2str(rho),',',num2str(theta),'),d2salphadthetadrho(',num2str(rho),',',num2str(theta),'),d2salphadrho2(',num2str(rho),',',num2str(theta),') = ',num2str(d2salphadtheta2),',',num2str(d2salphadthetadrho),',',num2str(d2salphadrho2)]);
1094         disp('----------------------------');
1095         disp(['gradrho(',num2str(rho),',',num2str(theta),'),dgradrhodrho(',num2str(rho),',',num2str(theta),'),dgradrhodtheta(',num2str(rho),',',num2str(theta),') = ',num2str(gradrho),',',num2str(dgradrhodrho),',',num2str(dgradrhodtheta)]);
1096         disp(['d2gradrhodtheta2(',num2str(rho),',',num2str(theta),'),d2gradrhodthetadrho(',num2str(rho),',',num2str(theta),'),d2gradrhodrho2(',num2str(rho),',',num2str(theta),') = ',num2str(d2gradrhodtheta2),',',num2str(d2gradrhodthetadrho),',',num2str(d2gradrhodrho2)]);
1097         disp('----------------------------');
1098         disp(['Upsilon(',num2str(rho),',',num2str(theta),'),dUpsilondrho(',num2str(rho),',',num2str(theta),'),dUpsilondtheta(',num2str(rho),',',num2str(theta),') = ',num2str(Upsilon),',',num2str(dUpsilondrho),',',num2str(dUpsilondtheta)]);
1099         disp(['d2Upsilondtheta2(',num2str(rho),',',num2str(theta),'),d2Upsilondthetadrho(',num2str(rho),',',num2str(theta),'),d2Upsilondrho2(',num2str(rho),',',num2str(theta),') = ',num2str(d2Upsilondtheta2),',',num2str(d2Upsilondthetadrho),',',num2str(d2Upsilondrho2)]);
1100    end
1101 end
1102