EAST_thomson_scattering_Tene_fit

 Calculates Te and ne profiles for C3PO/LUKE from Thomson scattering for EAST from raw data, using EQDSK equilibrium 

   INPUTS:

0001 function [equil_new,N_lav,rho_LUKE_ts_sr,rho_psi_ts_sr,T_ts_fit,N_ts_fit,T_ts_sr,dT_ts_sr,N_ts_sr,dN_ts_sr,T_rmse,N_rmse] = EAST_thomson_scattering_Tene_fit(equil,R_ts,Z_ts,T_ts,dT_ts,N_ts,dN_ts,R_hcn,display_mode,nefit_mode)
0002 %
0003 % Calculates Te and ne profiles for C3PO/LUKE from Thomson scattering for EAST from raw data, using EQDSK equilibrium
0004 %
0005 %   INPUTS:
0006 
0007 %       - equil: LUKE format magnetic equilibrium (struct)
0008 %       - R_ts: Major radius of the vertical chord of the diagnostic (m)
0009 %       - Z_ts: Vertical positions of the (ne,Te) raw measurements (m)
0010 %       - T_ts: Raw Te measurements (eV)
0011 %       - dT_ts: Standard deviation of raw Te measurements (eV)
0012 %       - N_ts: Raw ne measurements (m-3)
0013 %       - dN_ts: Standard deviation of raw ne measurements (m-3)
0014 %       - R_hcn: Major radius of the vertical chord of the HCN diagnostic that measures the line-averaged density (m)
0015 %       - display_mode: display data and fits (1: yes, 0: no)
0016 %           (default: 0)
0017 %       - nefit_mode: (1) bell-shaped curve guess function for the fit, (2) Fermi-like curve guess function for the fit
0018 %           (default: 2)
0019 %
0020 %   OUTPUTS:
0021 %
0022 %       - equil_new: LUKE magnetic equilibrium with new (ne,Te) profiles from Thomson scattering
0023 %       - N_lav: Line-averaged density (to be compared to HCN measurements) (m-3)
0024 %       - rho_LUKE_ts_sr: LUKE normalized radial coordinates (WARNING : not sqrt(poloidal_psin)) (a.u.)
0025 %       - rho_psi_ts_sr: rho as square root of normalized poloidal psi (a.u.)
0026 %       - T_ts_fit: Electron temperature profile from fit (keV)
0027 %       - N_ts_fit: Electron density profile from fit (m-3)
0028 %       - T_ts_sr: Electron temperature profile from Thomson scattering (keV)
0029 %       - dT_ts_sr: Standard deviation of electron temperature profile from Thomson scattering (keV)
0030 %       - N_ts_sr: Electron density profile from Thomson scattering (m-3)
0031 %       - dN_ts_sr: Standard deviation of electron density profile from Thomson scattering (m-3)
0032 %       - T_rmse: root mean square error for T
0033 %       - N_rmse: root mean square error for N
0034 %
0035 %   by Y. Peysson <yves.peysson@cea.fr> (CEA/IRFM)
0036 %
0037 if nargin == 8,
0038     display_mode = 0;
0039     nefit_mode = 2;
0040 end
0041 %
0042 if nargin == 9,
0043     nefit_mode = 2;
0044 end
0045 %
0046 if ~isfield(equil,'id')
0047     equil.id = '';
0048 end
0049 %
0050 rho_LUKE = equil.ptx(:,1)./equil.ptx(end,1);
0051 psin_LUKE = equil.psi_apRp./equil.psi_apRp(end);
0052 %
0053 equil_RZ = pt2rz_jd(equil);
0054 %
0055 T = @(x,a,b)a*exp(-x.^2/b^2);%Gaussian curve fit
0056 %
0057 if nefit_mode == 1    
0058     N = @(x,a,b)a*(1-x.*x).^b;%Bell shaped curvefit
0059 elseif nefit_mode == 2    
0060     N = @(x,a,b,c)a./(1+exp((x-b)/c));%Fermi function curve fit
0061 end
0062 %
0063 %
0064 % Sort data
0065 %
0066 [Z_ts_s,is] = sort(Z_ts,'descend');%Reorder vertical position
0067 T_ts_s = T_ts(is);
0068 dT_ts_s = dT_ts(is);
0069 %
0070 nonLinMdl_T_ts_s = fitnlm(Z_ts_s,T_ts_s,@(p,x)T(x,p(1),p(2)),[max(T_ts_s),1],'Weights',dT_ts_s);%fit raw Te data to cross-check the consistency between position of the magnetic axis from EFIT and raw Te-data
0071 Z_ts_s_ss = linspace(min(Z_ts_s),max(Z_ts_s),1000);%supergrid
0072 T_ts_s_ss_fit = T(Z_ts_s_ss,nonLinMdl_T_ts_s.Coefficients.Estimate('p1'),nonLinMdl_T_ts_s.Coefficients.Estimate('p2'));
0073 Zp_T_ts_s_ss_fit = Z_ts_s_ss(find(T_ts_s_ss_fit==max(T_ts_s_ss_fit)));
0074 %
0075 disp(['Zp (EFIT) = ',num2str(round(100*equil.Zp*100)/100),' (cm), Zp (Thomson Scattering) = ',num2str(round(100*Zp_T_ts_s_ss_fit*100)/100),' (cm)']);
0076 %
0077 N_ts_s = N_ts(is)/1e19;
0078 dN_ts_s = dN_ts(is)/1e19;
0079 N_ts_s(find(N_ts_s>100)) = NaN;% remove possible wrong points
0080 dN_ts_s(find(N_ts_s>100)) = NaN;% remove possible wrong points
0081 %
0082 X_ts_s = (R_ts - equil.Rp)*ones(size(Z_ts_s));
0083 Y_ts_s = Z_ts_s - equil.Zp;
0084 %
0085 rho_LUKE_ts_s = psi2rho_jd(equil,interp2(equil_RZ.x,equil_RZ.y,equil_RZ.xypsi_apRp'/equil.psi_apRp(end),X_ts_s,Y_ts_s));
0086 %
0087 Y_hcn = linspace(min(equil.pty(end,:)),max(equil.pty(end,:)),1000);
0088 grad_Y_hcn = gradient(Y_hcn);
0089 X_hcn = (R_hcn - equil.Rp)*ones(size(Y_hcn));
0090 %
0091 rho_hcn = psi2rho_jd(equil,interp2(equil_RZ.x,equil_RZ.y,equil_RZ.xypsi_apRp'/equil.psi_apRp(end),X_hcn,Y_hcn));
0092 %
0093 [rho_LUKE_ts_sr,isr] = sort(rho_LUKE_ts_s,'descend');
0094 T_ts_sr = T_ts_s(isr);
0095 dT_ts_sr = dT_ts_s(isr);
0096 N_ts_sr = N_ts_s(isr);
0097 dN_ts_sr = dN_ts_s(isr);
0098 %
0099 nonLinMdl_T_ts_sr = fitnlm(rho_LUKE_ts_sr,T_ts_sr,@(p,x)T(x,p(1),p(2)),[max(T_ts_sr),1],'Weights',dT_ts_sr);
0100 %
0101 if nefit_mode == 1       
0102     nonLinMdl_N_ts_sr = fitnlm(rho_LUKE_ts_sr,N_ts_sr,@(p,x)N(x,p(1),p(2)),[max(N_ts_sr),0.5],'Weights',dN_ts_sr);
0103 elseif nefit_mode == 2
0104     nonLinMdl_N_ts_sr = fitnlm(rho_LUKE_ts_sr,N_ts_sr,@(p,x)N(x,p(1),p(2),p(3)),[max(N_ts_sr),0.5,0.5],'Weights',dN_ts_sr);
0105 end
0106 %
0107 T_ts_fit = T(rho_LUKE,nonLinMdl_T_ts_sr.Coefficients.Estimate('p1'),nonLinMdl_T_ts_sr.Coefficients.Estimate('p2'));
0108 %
0109 if nefit_mode == 1
0110     N_ts_fit = N(rho_LUKE,nonLinMdl_N_ts_sr.Coefficients.Estimate('p1'),nonLinMdl_N_ts_sr.Coefficients.Estimate('p2'));
0111     N_hcn_fit = N(rho_hcn,nonLinMdl_N_ts_sr.Coefficients.Estimate('p1'),nonLinMdl_N_ts_sr.Coefficients.Estimate('p2'));
0112 elseif nefit_mode == 2
0113     N_ts_fit = N(rho_LUKE,nonLinMdl_N_ts_sr.Coefficients.Estimate('p1'),nonLinMdl_N_ts_sr.Coefficients.Estimate('p2'),nonLinMdl_N_ts_sr.Coefficients.Estimate('p3'));
0114     N_hcn_fit = N(rho_hcn,nonLinMdl_N_ts_sr.Coefficients.Estimate('p1'),nonLinMdl_N_ts_sr.Coefficients.Estimate('p2'),nonLinMdl_N_ts_sr.Coefficients.Estimate('p3'));
0115 end    
0116 %
0117 N_hcn_fit(isnan(N_hcn_fit)) = 0;%remove point outside the separatrix
0118 N_lav = sum(N_hcn_fit.*grad_Y_hcn)./sum(grad_Y_hcn);%line-averaged density to be compared to HCN measurements
0119 %
0120 if display_mode == 1  
0121     %
0122     figure,plot(Z_ts_s,T_ts_s,'ro',Z_ts_s_ss,T_ts_s_ss_fit,equil.Zp,T_ts_s_ss_fit(find(T_ts_s_ss_fit==max(T_ts_s_ss_fit))),'kx',Zp_T_ts_s_ss_fit,T_ts_s_ss_fit(find(T_ts_s_ss_fit==max(T_ts_s_ss_fit))),'k+');grid on
0123     hold on,errorbar(Z_ts_s,T_ts_s,dT_ts_s,'ro');
0124     axis([min(Z_ts) max(Z_ts) 0 max(T_ts_s)]);
0125     title([strrep(equil.id,'_','-'),', raw Te profile (Thomson scattering)']),xlabel(['Z@R=',num2str(R_ts),' (m)']),ylabel('T_e (eV)')
0126     legend('Exp.','Fit',['Z_p(EFIT) = ',num2str(round(100*equil.Zp*100)/100),' (cm)'],['Z_p(Fit) = ',num2str(round(100*Zp_T_ts_s_ss_fit*100)/100),' (cm)'])
0127     %
0128     figure,plot(Z_ts_s,N_ts_sr,'ro');grid on
0129     hold on,errorbar(Z_ts_s,N_ts_sr,dN_ts_sr,'ro');
0130     axis([min(Z_ts) max(Z_ts) 0 max(N_ts_sr)]);
0131     title([strrep(equil.id,'_','-'),', raw Ne profile (Thomson scattering)']),xlabel(['Z@R=',num2str(R_ts),' (m)']),ylabel('N_e (10^{+19} m^{-3})')
0132     %
0133     figure,plot(rho_LUKE_ts_sr,T_ts_sr,'ro',rho_LUKE,T_ts_fit);grid on
0134     hold on,errorbar(rho_LUKE_ts_sr,T_ts_sr,dT_ts_sr,'ro');
0135     axis([0 1 0 max(T_ts_sr)]);
0136     title([strrep(equil.id,'_','-'),', Te profile for LUKE']),xlabel('\rho_{LUKE}'),ylabel('T_e (eV)')
0137     legend('Exp.','Fit')
0138     %
0139     figure,plot(rho_LUKE_ts_sr,N_ts_sr,'ro',rho_LUKE,N_ts_fit);grid on
0140     hold on,errorbar(rho_LUKE_ts_sr,N_ts_sr,dN_ts_sr,'ro');
0141     axis([0 1 0 max(N_ts_sr)]);
0142     title([strrep(equil.id,'_','-'),', Ne profile for LUKE']),xlabel('\rho_{LUKE}'),ylabel('N_e (10^{+19} m^{-3})')
0143     legend('Exp.',['Fit, <n_e> = ',num2str(round(N_lav*1000)/1000),' 10^{+19} m^{-3}'],'Location','SouthWest')
0144     %
0145     figure,plot(sqrt(rho2psi_jd(equil,rho_LUKE_ts_sr)),T_ts_sr,'ro',sqrt(psin_LUKE),T_ts_fit);grid on
0146     hold on,errorbar(sqrt(rho2psi_jd(equil,rho_LUKE_ts_sr)),T_ts_sr,dT_ts_sr,'ro');
0147     axis([0 1 0 max(T_ts_sr)]);
0148     title([strrep(equil.id,'_','-'),', Te profile for LUKE']),xlabel('\rho = sqrt(\psi_{n})'),ylabel('T_e (eV)')
0149     legend('Exp.','Fit')
0150     %
0151     figure,plot(sqrt(rho2psi_jd(equil,rho_LUKE_ts_sr)),N_ts_sr,'ro',sqrt(psin_LUKE),N_ts_fit);grid on
0152     hold on,errorbar(sqrt(rho2psi_jd(equil,rho_LUKE_ts_sr)),N_ts_sr,dN_ts_sr,'ro');
0153     axis([0 1 0 max(N_ts_sr)]);
0154     title([strrep(equil.id,'_','-'),', Ne profile for LUKE']),xlabel('\rho = sqrt(\psi_{n})'),ylabel('N_e (10^{+19} m^{-3})')
0155     legend('Exp.',['Fit, <n_e> = ',num2str(round(N_lav*1000)/1000),' 10^{+19} m^{-3}'],'Location','SouthWest')
0156     %
0157     % for visual cross-check that 2-D interpolation works
0158     %
0159     %figure,plot(equil.ptx+equil.Rp,equil.pty+equil.Zp,'r',equil.ptx'+equil.Rp,equil.pty'+equil.Zp,'r',R_ts*ones(size(Z_ts)),Z_ts,'bo');
0160     %xx= equil_RZ.x'*ones(size(equil_RZ.y));
0161     %yy= equil_RZ.y'*ones(size(equil_RZ.x));
0162     %hold on,plot(xx+equil_RZ.Rp,yy'+equil_RZ.Zp,'k.');axis('equal');
0163     %contour(equil_RZ.x+equil_RZ.Rp,equil_RZ.y+equil_RZ.Zp,equil_RZ.xypsi_apRp'/equil.psi_apRp(end),100)
0164 %
0165 end
0166 %
0167 rho_psi_ts_sr = sqrt(rho2psi_jd(equil,rho_LUKE_ts_sr));%rho as square root of normalized poloidal psi
0168 %
0169 T_ts_sr = T_ts_sr/1000;%keV
0170 dT_ts_sr = dT_ts_sr/1000;%keV
0171 N_ts_sr = N_ts_sr*1e19;%m-3
0172 dN_ts_sr = dN_ts_sr*1e19;%m-3
0173 T_ts_fit = T_ts_fit/1000;%keV
0174 N_ts_fit = N_ts_fit*1e19;%m-3
0175 %
0176 T_rmse = nonLinMdl_T_ts_sr.RMSE/1000;
0177 N_rmse = nonLinMdl_N_ts_sr.RMSE*1e19;
0178 %
0179 %
0180 equil_new = equil;
0181 %
0182 equil_new.pTe = T_ts_fit;
0183 equil_new.pne = N_ts_fit;
0184 %
0185 
0186 
0187 
0188 
0189 
0190 
0191 
0192