load_externaldata_COMPASS

PURPOSE ^

SYNOPSIS ^

function external = load_externaldata_COMPASS(external_in,opt_gui)

DESCRIPTION ^

 This function loads data from COMPASS database

 INPUTS : 
   - if nargin == 0, return default external_in to be modified interactively in iload_externaldata_jd

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

0001 function external = load_externaldata_COMPASS(external_in,opt_gui)
0002     %
0003     % This function loads data from COMPASS database
0004     %
0005     % INPUTS :
0006     %   - if nargin == 0, return default external_in to be modified interactively in iload_externaldata_jd
0007 
0008     %   - external_in : input structure.
0009     %   - opt_gui     : option for GUI within function (interactive mode only)
0010     %
0011     %
0012     if nargin == 0,
0013         %
0014         % default values
0015         %
0016         external.shotnum = 0;
0017         external.shotime = 0;
0018         %
0019         % description & user halp
0020         %
0021         external.info.shotnum = 'Shot number';
0022         external.info.shotime = 'Shot time (s)';
0023         %
0024         % expert mode
0025         %
0026         external.expert.shotnum = 0;
0027         external.expert.shotime = 0;
0028         %
0029         % range of values
0030         %
0031         external.values.shotnum = [0;Inf];
0032         external.values.shotime = [0;Inf];
0033         %
0034         return
0035         %
0036     end
0037     %
0038     if nargin < 2,
0039         opt_gui = false;
0040     end
0041     %
0042     shotnum = external_in.shotnum;
0043     shotime = external_in.shotime;
0044     %
0045     % ******************************************************************
0046     %
0047     external = '';
0048     basestr = 'COMPASS';
0049     %
0050     if exist('cdb_client','file') == 2,%when we are on Abacus
0051         %
0052         data_out = compass2luke(shotnum);
0053         %
0054     else
0055         [filename,pathname] = igetfile_jd(opt_gui,'*.mat',['Please select the COMPASS data file for shot # ',shotnum],'',0);
0056         %
0057         if filename ~= 0,
0058             load([pathname,filename],'data_out');
0059         else
0060             disp(' ');disp('--> Read remotely from COMPASS database');disp(' ');
0061             [remote,profilestr] = load_remoteprofiles_jd('IPP_ABACUS',false,true,'direct',NaN,true);%read remotely from COMPASS database
0062             [err_data,~,data_out] = matremote_jd('compass2luke',{shotnum},remote,matstandardjob_jd('runmat'));        
0063             %
0064             if err_data ~= 0,
0065                 return
0066             end
0067         end
0068     end
0069     %
0070     if shotime == 0,% EFIT time grid reference, if shotime is not specified
0071         ishotime = round(length(data_out.efit.li.time_axis.data)/2);
0072         shotime = data_out.efit.li.time_axis.data(ishotime)/1e3;% ms -> s
0073         disp(' ');disp(['WARNING: The shot time has been enforced to ',num2str(shotime),' (s)']);disp(' ');
0074     else,%nearest value
0075         ishotime = interp1(data_out.efit.li.time_axis.data,[1:length(data_out.efit.li.time_axis.data)],shotime,'nearest')
0076         shotime = data_out.efit.li.time_axis.data(ishotime)/1e3;% ms -> s
0077         disp(' ');disp(['WARNING: The nearest shot time is ',num2str(shotime),' (s)']);disp(' ');
0078     end
0079     %
0080     % smoothing and interpolation (reference EFIT time grid)
0081     %
0082     %
0083     % EFIT data (see the definitions on the COMPASS website
0084     % https://www.tok.ipp.cas.cz/webcdb/data_sources/EFIT)
0085     %
0086     triangularity_lower_lcfs.t = data_out.efit.triangularity_lower_lcfs.time_axis.data(ishotime)/1e3;% ms -> s
0087     triangularity_lower_lcfs.data = data_out.efit.triangularity_lower_lcfs.data(ishotime);
0088     li.t = data_out.efit.li.time_axis.data(ishotime)/1e3;% ms -> s
0089     li.data = data_out.efit.li.data(ishotime);%internal inductance
0090     Ip.t = data_out.efit.Ip.time_axis.data(ishotime)/1e3;% ms -> s
0091     Ip.data = data_out.efit.Ip.data(ishotime);% A
0092     q.t = data_out.efit.q.time_axis.data(ishotime)/1e3;% ms -> s
0093     q.data = data_out.efit.q.data(ishotime,:);
0094     %
0095     psi_RZ.t = data_out.efit.psi_RZ.time_axis.data(ishotime)/1e3;% ms -> s
0096     psi_RZ.data = data_out.efit.psi_RZ.data(ishotime,:,:);% Wb/rad
0097     psi_RZ.R = data_out.efit.psi_RZ.axis1.data(ishotime,:);% m
0098     psi_RZ.Z = data_out.efit.psi_RZ.axis2.data(ishotime,:);% m
0099     psi_n.data = data_out.efit.psi_n.data;% Wb/rad
0100     psi_mag_axis.t = data_out.efit.psi_mag_axis.time_axis.data(ishotime)/1e3;% ms -> s
0101     psi_mag_axis.data = data_out.efit.psi_mag_axis.data(ishotime);
0102     psi_lcfs.t = data_out.efit.psi_lcfs.time_axis.data(ishotime)/1e3;% ms -> s
0103     psi_lcfs.data = data_out.efit.psi_lcfs.data(ishotime);
0104     %
0105     R_mag_axis.t = data_out.efit.R_mag_axis.time_axis.data(ishotime)/1e3;% ms -> s
0106     R_mag_axis.data = data_out.efit.R_mag_axis.data(ishotime,:,:);% m
0107     Z_mag_axis.t = data_out.efit.Z_mag_axis.time_axis.data(ishotime)/1e3;% ms -> s
0108     Z_mag_axis.data = data_out.efit.Z_mag_axis.data(ishotime,:,:);% m
0109     R_geom_axis.t = data_out.efit.R_geom_axis.time_axis.data(ishotime)/1e3;% ms -> s
0110     R_geom_axis.data = data_out.efit.R_geom_axis.data(ishotime,:,:);% m
0111     Z_geom_axis.t = data_out.efit.Z_geom_axis.time_axis.data(ishotime)/1e3;% ms -> s
0112     Z_geom_axis.data = data_out.efit.Z_geom_axis.data(ishotime,:,:);% m
0113     %
0114     B_vac_R_geom.t = data_out.efit.B_vac_R_geom.time_axis.data(ishotime)/1e3;% ms -> s
0115     B_vac_R_geom.data = data_out.efit.B_vac_R_geom.data(ishotime,:,:);% m
0116     RBphi.t = data_out.efit.RBphi.time_axis.data(ishotime)/1e3;% ms -> s
0117     RBphi.data = data_out.efit.RBphi.data(ishotime,:,:);% m
0118     FFprime.t = data_out.efit.FFprime.time_axis.data(ishotime)/1e3;% ms -> s
0119     FFprime.data = data_out.efit.FFprime.data(ishotime,:,:);% m
0120     pprime_static.t = data_out.efit.pprime_static.time_axis.data(ishotime)/1e3;% ms -> s
0121     pprime_static.data = data_out.efit.pprime_static.data(ishotime,:,:);% m
0122     p_stat.t = data_out.efit.p_stat.time_axis.data(ishotime)/1e3;% ms -> s
0123     p_stat.data = data_out.efit.p_stat.data(ishotime,:,:);% m
0124     %
0125     if isfield(data_out,'ne') && isfield(data_out,'te'),%Thomson scattering
0126         pne.t = shotime;
0127         pne.x = data_out.ne.axis1.data;%Vertical TS measuring points (m)
0128         pne.data = interp1(data_out.ne.time_axis.data/1e3,data_out.ne.data,shotime,'linear');% m^-3
0129         pne.data = compass_ts_profile_fit(pne,psi_RZ,psi_n.data,psi_lcfs.data,psi_mag_axis.data);% m^-3
0130         %
0131         pTe.t = shotime;
0132         pTe.x = data_out.te.axis1.data;%Vertical TS measuring points (m)
0133         pTe.data = interp1(data_out.te.time_axis.data/1e3,data_out.te.data/1e3,shotime,'linear');% keV
0134         pTe.data = compass_ts_profile_fit(pTe,psi_RZ,psi_n.data,psi_lcfs.data,psi_mag_axis.data);% keV
0135     else
0136         pne.t = [];
0137         pne.data = [];
0138         pTe.t = [];
0139         pTe.data = [];
0140     end
0141     %
0142     % Magnetic data
0143     %
0144     U_loop.t = shotime;
0145     %U_loop.data = interp1(data_out.U_loop.time_axis.data/1e3,sgolayfilt(data_out.U_loop.data, 1,min(2*round(length(data_out.U_loop.time_axis.data)/400)-1,10001)),shotime,'linear');% V
0146     t_resample = 100;
0147     U_loop.data = interp1(data_out.U_loop.time_axis.data(t_resample/2:t_resample:end-t_resample/2)/1e3,compass_resample_data(data_out.U_loop.data,t_resample),shotime,'linear');% V
0148     %close all,plot(data_out.U_loop.time_axis.data/1e3,data_out.U_loop.data,'b',data_out.U_loop.time_axis.data/1e3,sgolayfilt(data_out.U_loop.data, 1,min(2*round(length(data_out.U_loop.time_axis.data)/400)-1,10001)),'r')
0149     I_plasma.t = shotime;
0150     %I_plasma.data = interp1(data_out.I_plasma.time_axis.data/1e3,sgolayfilt(data_out.I_plasma.data, 1,min(2*round(length(data_out.I_plasma.time_axis.data)/400)-1,10001)),shotime,'linear');% A
0151     I_plasma.data = interp1(data_out.I_plasma.time_axis.data(t_resample/2:t_resample:end-t_resample/2)/1e3,compass_resample_data(data_out.I_plasma.data,t_resample),shotime,'linear');% A
0152     %close all,plot(data_out.I_plasma.time_axis.data/1e3,data_out.I_plasma.data,'b',data_out.I_plasma.time_axis.data/1e3,sgolayfilt(data_out.I_plasma.data, 1,min(2*round(length(data_out.I_plasma.time_axis.data)/400)-1,10001)),'r')
0153     %
0154     % return structure external
0155     %
0156     external.basestr = basestr;
0157     external.shotnum = num2str(shotnum);
0158     external.shotime = shotime;
0159     %
0160     external.ohm.vloop = U_loop.data;
0161     %
0162     external.equil.shotnum = external.shotnum;
0163     external.equil.shotime = external.shotime;
0164     %
0165     external.equil.magnetic.I_plasma = I_plasma.data;
0166     %
0167     external.equil.magnetic.efit.triangularity_lower_lcfs = triangularity_lower_lcfs.data;
0168     external.equil.magnetic.efit.li = li.data;
0169     external.equil.magnetic.efit.Ip = Ip.data;
0170     external.equil.magnetic.efit.q = q.data;
0171     %
0172     external.equil.magnetic.efit.psi_RZ = reshape(psi_RZ.data,[length(psi_RZ.R),length(psi_RZ.Z)]);% Wb/rad
0173     external.equil.magnetic.efit.psi_n = psi_n.data;
0174     external.equil.magnetic.efit.psi_mag_axis = psi_mag_axis.data;
0175     external.equil.magnetic.efit.psi_lcfs = psi_lcfs.data;
0176     %
0177     external.equil.magnetic.efit.R = psi_RZ.R;% m
0178     external.equil.magnetic.efit.Z = psi_RZ.Z;% m
0179     external.equil.magnetic.efit.R_mag_axis = R_mag_axis.data;% m
0180     external.equil.magnetic.efit.Z_mag_axis = Z_mag_axis.data;% m
0181     external.equil.magnetic.efit.R_geom_axis = R_geom_axis.data;% m
0182     external.equil.magnetic.efit.Z_geom_axis = Z_geom_axis.data;% m
0183     %
0184     external.equil.magnetic.efit.B_vac_R_geom = B_vac_R_geom.data;% m
0185     external.equil.magnetic.efit.RBphi = RBphi.data;% m
0186     external.equil.magnetic.efit.FFprime = FFprime.data;% m
0187     external.equil.magnetic.efit.pprime_static = pprime_static.data;% m
0188     external.equil.magnetic.efit.p_stat = p_stat.data;% m
0189     %
0190     external.equil.prof.pne = pne.data;
0191     external.equil.prof.pTe = pTe.data;
0192     %
0193     external.equil.prof.pspsin = psi_n.data;
0194     external.equil.prof.pne = pne.data;
0195     external.equil.prof.pTe = pTe.data;
0196     external.equil.prof.pTi = pTe.data/2;%Ti = Te/2
0197     external.equil.prof.Zeff = 2;%Zeff = 2
0198     external.equil.prof.fi = [0,1,0];%Deuterium plasma only
0199     external.equil.prof.zZi = [1,1,1,6];%carbon main impurity
0200     external.equil.prof.zmi = [1,2,3,12];%carbon main impurity
0201     %
0202     external.id = ['COMPASS_',external.shotnum,'_',num2str(shotime,'%6.4f')];
0203     %
0204     if isnan_jd(external,0),
0205         disp('Warning : NaNs were found in the structure ''external'' : ')
0206     end
0207     %
0208 end
0209 
0210 function [data_out] = compass_ts_profile_fit(data_in,psi_RZ,psi_n,psi_lcfs,psi_mag_axis)
0211     %
0212     % Thomson scattering profile smoothing (remove spurious points)
0213     %
0214     x = data_in.x;
0215     ts = data_in.data;
0216     %
0217     %psi1  = psi_lcfs;
0218     %psi0  = psi_mag_axis;
0219     %
0220     psi_ts = interp2(psi_RZ.R,psi_RZ.Z,squeeze(psi_RZ.data).',0.56,x,'spline');%psi value where TS measures
0221     psin_ts = (psi_ts - psi_mag_axis)/(psi_lcfs - psi_mag_axis);%normalized psi where TS measures
0222     %
0223     ipoints = find(x < 0.18);%0.33
0224     %
0225     ts_ptot = ts(ipoints);
0226     %
0227     figure,subplot(2,1,1)
0228     title('Thomson scattering density profile')     
0229     hold on     
0230     plot(x(ipoints),ts(ipoints) ,'rx');
0231     plot(x,ts,'b');
0232     %     ylim([0 4e19]);
0233     ylabel('Profile values');
0234     xlabel('Z coordinate [m]');
0235     hold off
0236     %
0237     ts_psin = psin_ts(ipoints); 
0238     %
0239     sfactor = 0.999;     
0240     ptot_fit = fit(psin_ts(ipoints), ts_ptot.', 'smoothingspline', 'SmoothingParam', sfactor);
0241     data_out = ptot_fit(psi_n);
0242     %
0243     subplot(2,1,2) % first subplot
0244     hold on
0245     plot(psin_ts(ipoints), ts_ptot.', 'rx')
0246     plot(psi_n, data_out, 'b')
0247     %     ylim([0 4e19])
0248     ylabel('Profile values')
0249     xlabel('Normalized \psi')
0250     hold off
0251     %
0252 end
0253 
0254 function res = compass_resample_data(data, ww)
0255 
0256     % fast resample data using a running average filter
0257     %
0258     % data: vector of data
0259     % ww: window width
0260 
0261     % enforce column vector
0262     % TODO adjust automatically for row/column vectors
0263     data = reshape(data, [], 1);
0264     % cut the last data points
0265     % TODO properly treat the last data points
0266     n_mod = mod(length(data), ww);
0267     data = data(1:end-n_mod);
0268     % reshape the data and calculate the mean
0269     data = reshape(data, ww, []);
0270     res = mean(data, 1);
0271 
0272 end
0273
Community support and wiki are available on Redmine. Last update: 18-Apr-2019.