testf0struct_yp

 Test the interpolated distribution function and PSI =  B/B0 where B0 is
 the minimum of B on a magnetic flux surface

 INPUT

0001 function [PSI_fit1,Xf0_fit] = testf0struct_yp(dke_out,momentumDKE,equilDKE,mksa,f0struct,rho,theta,pn,mhu)
0002 %
0003 % Test the interpolated distribution function and PSI =  B/B0 where B0 is
0004 % the minimum of B on a magnetic flux surface
0005 %
0006 % INPUT
0007      
0008 %    - pn: momentum normalized to pth_ref (betath_ref) [1,1]
0009 %      (default = 4)
0010 %    - mhu: pitch-angle (-1 -> 1) [1,1]
0011 %      (default = 0.5)
0012 %    - rho: normalized minor radius [1,1]
0013 %      (default = 0.25)
0014 %    - theta: poloidal angle (rad) [1,1]
0015 %      (default = 0.2)
0016 %
0017 % OUTPUT:
0018 %    - PSI: B/B0 where B0 is the minimum of B on a magnetic flux surface [1,1]
0019 %    - Xf0: distribution function at B0(rho) [npn,nmhu]
0020 %
0021 % By Yves Peysson (CEA-DRFC, yves.peysson@cea.fr) and Joan Decker (CEA-DRFC, joan.decker@cea.fr)
0022 %
0023 if nargin < 5,
0024     error('Wrong number of input arguments in testf0struct_yp.m');
0025 end
0026 %
0027 if nargin < 6,%Output for comparison with C raytracing_yp.c
0028    rho = 0.25;
0029    theta = 0.2;
0030    pn = 4;
0031    mhu = 0.5;
0032 end
0033 %
0034 if nargin < 7,
0035    theta = 0.2;
0036    pn = 4;
0037    mhu = 0.5;
0038 end
0039 %
0040 if nargin < 8,
0041    pn = 4;
0042    mhu = 0.5;
0043 end
0044 %
0045 if nargin < 9,
0046    mhu = 0.5;
0047 end
0048 %
0049 ipn = max(find(f0struct.pn <= pn));
0050 jmhu = max(find(f0struct.mhu <= mhu));
0051 pnd = f0struct.pn(ipn);
0052 mhud = f0struct.mhu(jmhu);
0053 %
0054 B = sqrt(equilDKE.XBx.*equilDKE.XBx + equilDKE.XBy.*equilDKE.XBy + equilDKE.XBphi.*equilDKE.XBphi);%Module of the total magnetic field
0055 PSI = B./(B(:,1)*ones(1,length(equilDKE.mtheta)));%Module of the total magnetic field at theta = 0 (where it is minimum)
0056 %
0057 if nargin < 6, 
0058     close all;
0059     %
0060     ir = 4;%resampling factor
0061     %
0062     rho_fit = sort(resample(equilDKE.xrho,ir,1)');
0063     for ir = 1:length(equilDKE.xrho),
0064         index_rho(ir) = find(rho_fit == equilDKE.xrho(ir));
0065     end
0066     theta_fit = equilDKE.mtheta';
0067     %
0068     PSI_a0_fit = ppval(f0struct.PSI_fit.pp_a0,rho_fit);
0069     PSI_an_fit = ppval(f0struct.PSI_fit.pp_an,rho_fit);
0070     PSI_bn_fit = ppval(f0struct.PSI_fit.pp_bn,rho_fit);
0071     PSI_da0drho_fit = ppval(f0struct.PSI_fit.pp_da0drho,rho_fit);
0072     PSI_dandrho_fit = ppval(f0struct.PSI_fit.pp_dandrho,rho_fit);
0073     PSI_dbndrho_fit = ppval(f0struct.PSI_fit.pp_dbndrho,rho_fit);
0074     %
0075     PSI_a0_fit1 = ppval_yp(f0struct.PSI_fit.pp_a0,rho_fit);
0076     PSI_an_fit1 = ppval_yp(f0struct.PSI_fit.pp_an,rho_fit);
0077     PSI_bn_fit1 = ppval_yp(f0struct.PSI_fit.pp_bn,rho_fit);
0078     PSI_da0drho_fit1 = ppval_yp(f0struct.PSI_fit.pp_da0drho,rho_fit);
0079     PSI_dandrho_fit1 = ppval_yp(f0struct.PSI_fit.pp_dandrho,rho_fit);
0080     PSI_dbndrho_fit1 = ppval_yp(f0struct.PSI_fit.pp_dbndrho,rho_fit);
0081     %
0082     % Build interpolated quantities
0083     %
0084     PSI_fit = zeros(length(rho_fit),length(theta_fit));
0085     dPSIdtheta_fit = zeros(length(rho_fit),length(theta_fit));
0086     d2PSIdtheta2_fit = zeros(length(rho_fit),length(theta_fit));
0087     dPSIdrho_fit = zeros(length(rho_fit),length(theta_fit));
0088     d2PSIdthetadrho_fit = zeros(length(rho_fit),length(theta_fit));
0089     %
0090     PSI_fit1 = zeros(length(rho_fit),length(theta_fit));
0091     dPSIdtheta_fit1 = zeros(length(rho_fit),length(theta_fit));
0092     d2PSIdtheta2_fit1 = zeros(length(rho_fit),length(theta_fit));
0093     dPSIdrho_fit1 = zeros(length(rho_fit),length(theta_fit));
0094     d2PSIdthetadrho_fit1 = zeros(length(rho_fit),length(theta_fit));
0095     %
0096     for ii = 1:length(theta_fit), 
0097         [PSI_fit(:,ii)] = calcval_yp(f0struct.PSI_fit,theta_fit(ii),PSI_a0_fit(:),PSI_an_fit,PSI_bn_fit,PSI_da0drho_fit(:),PSI_dandrho_fit,PSI_dbndrho_fit);
0098         [PSI_fit1(:,ii)] = calcval_yp(f0struct.PSI_fit,theta_fit(ii),PSI_a0_fit1(:),PSI_an_fit1,PSI_bn_fit1,PSI_da0drho_fit1(:),PSI_dandrho_fit1,PSI_dbndrho_fit1);
0099     end
0100     %
0101     figure('Name','PSI'),
0102     [ax] = graph1D_jd(theta_fit,PSI_fit(index_rho,:),0,0,'','','',NaN,NaN,NaN,'-','none','b',2);
0103     [ax] = graph1D_jd(theta_fit,PSI_fit1(index_rho,:),0,0,'','','',NaN,NaN,NaN,'-','none','c',1);    
0104     [ax] = graph1D_jd(equilDKE.mtheta,PSI,0,0,'\theta (rd)','\Psi','\Psi(\rho,\theta) = B(\rho,\theta)/B_{0}(\rho)',NaN,[0,2*pi],NaN,'--','none','r',1);
0105     legend('Matlab interp. (blue)','LUKE interp. (cyan)','LUKE (red)')
0106     %
0107     Xf0_fit = reshape(ppval(f0struct.XXf0_fit.pp_f,rho_fit)',length(f0struct.pn),length(f0struct.mhu),length(rho_fit));
0108     ipn = max(find(f0struct.pn <= pn));
0109     %
0110     ne_fit = interp1(equilDKE.xrho,equilDKE.xne,rho_fit)/mksa.ne_ref;%interpolation fo the electron density profile
0111     %
0112     Xpn2_fit = (f0struct.pn(:).^2)*ones(1,length(f0struct.mhu));
0113     for ir = 1:length(rho_fit)
0114         normXf0_fit(ir) = 2*pi*integral_dke_jd(momentumDKE.dmhu,integral_dke_jd(momentumDKE.dpn,Xf0_fit(:,:,ir).*Xpn2_fit,1),1)';
0115     end
0116     %
0117     figure('Name','f0'),
0118     [ax] = graph1D_jd(rho_fit,squeeze(Xf0_fit(ipn,jmhu,:)),0,0,'','','',NaN,NaN,NaN,'-','none','b',2);
0119     [ax] = graph1D_jd(equilDKE.xrho,squeeze(dke_out.XXf0(ipn,jmhu,:)),0,0,'\rho','f_{0}',['@ p_{n} = ',num2str(pnd),', \xi_{0} = ',num2str(mhud)],NaN,[0,1],NaN,'none','+','r',2);
0120     legend('Fit','LUKE')
0121     %
0122     figure('Name','Norm. f0'),
0123     [ax] = graph1D_jd(rho_fit,normXf0_fit,0,0,'','','',NaN,NaN,NaN,'-','none','b',2);
0124     [ax] = graph1D_jd(rho_fit,ne_fit,0,0,'\rho','','',NaN,NaN,NaN,'-','none','r',1);
0125     [ax] = graph1D_jd(equilDKE.xrho,equilDKE.xne/mksa.ne_ref,0,0,'\rho','Normalized electron density','',NaN,[0,1],NaN,'none','+','r',2);
0126     legend('Fit of f_{0}','Fit of n_{e}','n_{e}')
0127 else
0128     %
0129     PSI_a0_fit1 = ppval_yp(f0struct.PSI_fit.pp_a0,rho);
0130     PSI_an_fit1 = ppval_yp(f0struct.PSI_fit.pp_an,rho);
0131     PSI_bn_fit1 = ppval_yp(f0struct.PSI_fit.pp_bn,rho);
0132     PSI_da0drho_fit1 = ppval_yp(f0struct.PSI_fit.pp_da0drho,rho);
0133     PSI_dandrho_fit1 = ppval_yp(f0struct.PSI_fit.pp_dandrho,rho);
0134     PSI_dbndrho_fit1 = ppval_yp(f0struct.PSI_fit.pp_dbndrho,rho);
0135     %
0136     [PSI_fit1] = calcval_yp(f0struct.PSI_fit,theta,PSI_a0_fit1(:),PSI_an_fit1,PSI_bn_fit1,PSI_da0drho_fit1(:),PSI_dandrho_fit1,PSI_dbndrho_fit1);
0137     %
0138     Xf0_fit = reshape(ppval(f0struct.XXf0_fit.pp_f,rho)',length(f0struct.pn),length(f0struct.mhu),length(rho));
0139     ne_fit = interp1(equilDKE.xrho,equilDKE.xne,rho)/mksa.ne_ref;%interpolation fo the electron density profile
0140     %
0141     Xpn2_fit = (f0struct.pn(:).^2)*ones(1,length(f0struct.mhu));
0142     normXf0_fit = 2*pi*integral_dke_jd(momentumDKE.dmhu,integral_dke_jd(momentumDKE.dpn,Xf0_fit.*Xpn2_fit,1),1)';
0143     %
0144     disp('---------------------------------');
0145     disp(['XXf0_fit.pp_f.coefs(1) = ',num2str(f0struct.XXf0_fit.pp_f.coefs(1))]);
0146     disp(['XXf0_fit.pp_dfdrho.coefs(1) = ',num2str(f0struct.XXf0_fit.pp_dfdrho.coefs(1))]);
0147     disp(['XXf0_fit.pp_d2fdrho2.coefs(1) = ',num2str(f0struct.XXf0_fit.pp_d2fdrho2.coefs(1))]);
0148     disp('--------------');        
0149     disp(['Xf0_fit(',num2str(pnd),',',num2str(mhud),',',num2str(rho),') = ',num2str(squeeze(Xf0_fit(ipn,jmhu)))]);
0150     disp(['ne/ne_ref(',num2str(rho),') = ',num2str(ne_fit)]);
0151     disp(['norm_Xf0_fit(',num2str(rho),') = ',num2str(normXf0_fit)]);
0152     disp('---------------------------------');      
0153     disp(['PSI_fit.n = ',num2str(f0struct.PSI_fit.n(end))]);
0154     disp(['PSI_fit.pp_bn.coefs(1) = ',num2str(f0struct.PSI_fit.pp_bn.coefs(1))]);
0155     disp(['PSI_fit.pp_dbndrho.coefs(1) = ',num2str(f0struct.PSI_fit.pp_dbndrho.coefs(1))]);
0156     disp('--------------');
0157     disp(['PSI(',num2str(rho),',',num2str(theta),') = ',num2str(PSI_fit1)]);
0158     disp('---------------------------------');   
0159 end
0160