fproc_fluctn_jd

0001 function Zfluct_out = fproc_fluctn_jd(filename,fluct_display,Zfluct_ref,leg_ref,opt)
0002 %
0003 if nargin < 5,
0004     opt = struct;
0005 end
0006 %
0007 if nargin < 3,
0008     Zfluct_ref = '';
0009 end
0010 %
0011 if nargin < 2 || isempty(fluct_display),
0012     fluct_display.mode = 0;
0013     fluct_display.evol = 0;
0014     fluct_display.stat = 0;
0015 end
0016 %
0017 if nargin < 4,
0018     nref = length(Zfluct_ref);
0019     if nref == 1,
0020         leg_ref = {'ref.'};
0021     else
0022         leg_ref = cell(1,nref);
0023         for iref=1:nref,
0024             leg_ref{iref} = ['ref. ',num2str(iref)];
0025         end
0026     end
0027 end
0028 %
0029 if isstruct(filename),% iluke data
0030     %
0031     lukedata = load(filename.filename);
0032     %
0033     ishot = find(lukedata.shotimes == filename.shotime);
0034     output = lukedata.outputs{ishot};
0035     data.equil = lukedata.equils{ishot};
0036     %
0037     data.Ztime = output.dke_out.FPtime;
0038     data.Zitn = length(data.Ztime);
0039     data.Zcurr = output.Zcurr;
0040     data.ZP0 = output.ZP0;
0041     data.XXf0_t = output.dke_out.XXf0;
0042     data.XXfM = output.dke_out.XXfM;
0043     data.simul = output.dke_out.simul;
0044     data.equil = output.dke_out.equil;
0045     equilDKE = output.equilDKE;
0046     mksa = output.mksa;
0047     radialDKE = output.radialDKE;
0048     momentumDKE = output.momentumDKE;
0049     %
0050     clear lukedata ishot output 
0051     %
0052 else
0053     %
0054     data = load(filename,'Ztime','Zitn','Zcurr','ZP0','dke_out','equil','XXf0_t','simul','loadstructs','calctime','lukestructs');
0055     %
0056     if isfield(data,'loadstructs') && isfield(data.loadstructs,'equilibrium_jd'),% new run_lukert, backup file, with fluct
0057         equilDKE = data.loadstructs.equilibrium_jd.equilDKE;
0058         mksa = data.loadstructs.mksacoefbuilder_dke_yp.mksa;
0059         radialDKE = data.loadstructs.radialgrid_dke_jd.radialDKE;
0060         momentumDKE = data.loadstructs.momentumgrid_dke_jd.momentumDKE;
0061         data.XXfM = data.loadstructs.coll_dke_jd.XXfM;
0062     elseif isfield(data,'loadstructs'),% old run_lukert, backup file, with fluct
0063         equilDKE = data.loadstructs.equilDKE;
0064         mksa = data.loadstructs.mksa;
0065         radialDKE = data.loadstructs.radialDKE;
0066         momentumDKE = data.loadstructs.momentumDKE;
0067         data.XXfM = data.loadstructs.XXfM;
0068     elseif isfield(data,'lukestructs'),% new run_lukert, LUKE_RESULTS, no fluct
0069         data.Ztime = data.lukestructs.output.dke_out.FPtime;
0070         data.Zitn = length(data.Ztime);
0071         data.Zcurr = data.lukestructs.output.Zcurr;
0072         data.ZP0 = data.lukestructs.output.ZP0;
0073         data.XXf0_t = data.lukestructs.output.dke_out.XXf0;
0074         data.XXfM = data.lukestructs.output.dke_out.XXfM;
0075         data.simul = data.lukestructs.output.dke_out.simul;
0076         data.equil = data.lukestructs.output.dke_out.equil;
0077         equilDKE = data.lukestructs.output.equilDKE;
0078         mksa = data.lukestructs.output.mksa;
0079         radialDKE = data.lukestructs.output.radialDKE;
0080         momentumDKE = data.lukestructs.output.momentumDKE;
0081     else% old run_lukert, LUKE_RESULT, no fluct
0082         data = conc_struct_jd(data,load(filename,'equilDKE','mksa','radialDKE','momentumDKE'));
0083         data.XXf0_t = data.dke_out.XXf0;
0084         data.XXfM = data.dke_out.XXfM;
0085         data.simul = data.dke_out.simul;
0086         data.equil = data.dke_out.equil;
0087         equilDKE = data.equilDKE;
0088         mksa = data.mksa;
0089         radialDKE = data.radialDKE;
0090         momentumDKE = data.momentumDKE;
0091     end
0092 end
0093 %
0094 xrho = equilDKE.xrho.';%geometrical rho (xrhoG in proc_luke_jd.m)
0095 xdrho = diff(radialDKE.xrho_S_dke).';
0096 %
0097 xrhoT = equilDKE.xrhoT.'; 
0098 xrhoT_S = [0,interp1(xrho,xrhoT,radialDKE.xrho_S_dke(2:end-1),'spline'),1];
0099 xdrhoT = diff(xrhoT_S.');
0100 %
0101 fitrhoT = fit(xrho,xrhoT,'cubicinterp');
0102 %
0103 xdA = equilDKE.xdA_dke.';
0104 %
0105 if isfield(equilDKE,'xdV_dke'),
0106     xdV = equilDKE.xdV_dke.'; 
0107 else
0108     xdV = equilDKE.xdV_2piRp_dke.'*2*pi*equilDKE.Rp;
0109 end
0110 %
0111 if isfield(data,'Ztime'),% new run_lukert
0112     Ztime = data.Ztime;%cumsum(data.Ztime);%
0113 elseif isfield(data.Zcurr,'xtn'),% old run_lukert, with fluct
0114     Ztime = data.Zcurr.xtn;
0115 else% old run_lukert, no fluct
0116     data = conc_struct_jd(data,load(filename,'dke_out'));
0117     Ztime = data.dke_out.FPtime;
0118 end
0119 %
0120 nt = length(Ztime);
0121 nr = length(xrho);
0122 %
0123 if isfield(data,'Zitn'),% new run_lukert
0124     Zitn = data.Zitn;
0125 elseif isfield(data.Zcurr,'xtn'),% old run_lukert, with fluct
0126     Zitn = 1:nt;
0127 else% old run_lukert, no fluct
0128     Zitn = nt;
0129 end
0130 %
0131 nitn = length(Zitn);
0132 %
0133 if isfield(data,'Ztime') && any(diff(Ztime) <= 0),% keep until corrected in run_lukert
0134     Zcumtime = [0,cumsum(Ztime(Zitn))];
0135     Ztime = Ztime + Zcumtime(ceil((1:nt)/Zitn(1)));
0136 end
0137 %
0138 calctime.waves = NaN(1,nitn);
0139 calctime.luke = NaN(1,nitn);
0140 %
0141 if isfield(data,'calctime'),
0142     for itn = 1:nitn,
0143         calctime.waves(itn) = data.calctime(itn).waves;
0144         calctime.luke(itn) = data.calctime(itn).luke;
0145     end
0146 end
0147 %
0148 if length(data.Zcurr) == nt,% new run_lukert
0149     J = NaN(nt,nr);
0150     P = NaN(nt,nr);
0151     %
0152     for it = 1:nt,
0153         %
0154         J(it,:) = data.Zcurr(it).x_0_fsav*mksa.j_ref;
0155         P(it,:) = data.ZP0(it).x_rf_fsav*mksa.P_ref;
0156     end
0157 elseif isfield(data.Zcurr,'xtn'),% old run_lukert, with fluct
0158     J = data.Zcurr.t_x_0_fsav*mksa.j_ref;
0159     P = data.ZP0.t_x_rf_fsav*mksa.P_ref;
0160 end
0161 %
0162 I_tot = J*xdA;
0163 P_rf = P*xdV;
0164 %
0165 if nitn > 1 && fluct_display.evol,
0166     plot_fluct_jd(Zitn,xrho,momentumDKE.pn,J,P,data.XXfM,data.XXf0_t,data.simul,data.equil.pTe(1),opt)
0167 end
0168 %
0169 %
0170 %
0171 fwhm = 2*sqrt(2*log(2));% ratio of fwhm to STD
0172 %
0173 Zmask = Zitn;%
0174 %
0175 if isfield(opt,'racc'),
0176     racc = opt.racc;
0177 else
0178     racc.P = NaN;
0179     racc.I = NaN;
0180 end
0181 %
0182 if isfield(opt,'itnmin'),
0183     itnmin = opt.itnmin;
0184 else
0185     itnmin = 1;
0186 end
0187 %
0188 [~,~,fluctres] = fluctanalysis_yp(Ztime(Zmask),J(Zmask,:),P(Zmask,:),I_tot(Zmask),P_rf(Zmask),racc,itnmin);
0189 %
0190 sitn = (itnmin - 1) + fluctres.P.iPtot(fluctres.I.iItot);%select values after rejection of abnormal power and current ramp-up points
0191 sitnPr = [1:itnmin - 1, (itnmin - 1) + fluctres.P.iPtotr];
0192 sitnIr = (itnmin - 1) + fluctres.P.iPtot(fluctres.I.iItotr);
0193 %
0194 stn = Ztime(Zmask(sitn));
0195 %
0196 ntn = length(stn);
0197 %
0198 sJ = J(Zmask(sitn),:);% time-dependent current profile
0199 sP = P(Zmask(sitn),:);% time-dependent power profile
0200 sI_tot = I_tot(Zmask(sitn));% time-dependent current
0201 sP_rf = P_rf(Zmask(sitn));% time-dependent powert
0202 %
0203 srhoJ = sJ*(xrho.*xdA)./sI_tot;% time-dependent current peak location
0204 srhoP = sP*(xrho.*xdV)./sP_rf;% time-dependent power peak location
0205 %
0206 srhoJT = sJ*(xrhoT.*xdA)./sI_tot;% time-dependent current peak location rhoT grid
0207 srhoPT = sP*(xrhoT.*xdV)./sP_rf;% time-dependent power peak location rhoT grid
0208 %
0209 srho = repmat(xrho.',[ntn,1]);
0210 srhoT = repmat(xrhoT.',[ntn,1]);%rhoT grid
0211 sxrhoJ = repmat(srhoJ,[1,nr]);
0212 sxrhoJT = repmat(srhoJT,[1,nr]);%rhoT grid
0213 sxrhoP = repmat(srhoP,[1,nr]);
0214 sxrhoPT = repmat(srhoPT,[1,nr]);%rhoT grid
0215 %
0216 %
0217 sdrhoJ = fwhm*sqrt((sJ.*((srho-sxrhoJ).^2)*xdA)./sI_tot);%time-dependent current FWHM
0218 sdrhoP = fwhm*sqrt((sP.*((srho-sxrhoP).^2)*xdV)./sP_rf);%time-dependent power FWHM
0219 %
0220 sdrhoJT = fwhm*sqrt((sJ.*((srhoT-sxrhoJT).^2)*xdA)./sI_tot);%time-dependent current FWHM rhoT grid
0221 sdrhoPT = fwhm*sqrt((sP.*((srhoT-sxrhoPT).^2)*xdV)./sP_rf);%time-dependent power FWHM rhoT grid
0222 %
0223 sI_totm = NaN(1,ntn);
0224 ssI_tot =  NaN(1,ntn);
0225 srI_totm =  NaN(1,ntn);
0226 %
0227 sP_rfm =  NaN(1,ntn);
0228 ssP_rf =  NaN(1,ntn);
0229 srP_rfm =  NaN(1,ntn);
0230 %
0231 srhoJm =  NaN(1,ntn);
0232 ssrhoJ =  NaN(1,ntn);
0233 srrhoJm =  NaN(1,ntn);
0234 %
0235 srhoJmT =  NaN(1,ntn);
0236 ssrhoJT =  NaN(1,ntn);
0237 srrhoJmT =  NaN(1,ntn);
0238 %
0239 srhoPm =  NaN(1,ntn);
0240 ssrhoP =  NaN(1,ntn);
0241 srrhoPm =  NaN(1,ntn);
0242 %
0243 srhoPmT =  NaN(1,ntn);
0244 ssrhoPT =  NaN(1,ntn);
0245 srrhoPmT =  NaN(1,ntn);
0246 %
0247 sdrhoJm =  NaN(1,ntn);
0248 ssdrhoJ =  NaN(1,ntn);
0249 srdrhoJm =  NaN(1,ntn);
0250 %
0251 sdrhoJmT =  NaN(1,ntn);
0252 ssdrhoJT =  NaN(1,ntn);
0253 srdrhoJmT =  NaN(1,ntn);
0254 %
0255 sdrhoPm =  NaN(1,ntn);
0256 ssdrhoP =  NaN(1,ntn);
0257 srdrhoPm =  NaN(1,ntn);
0258 %
0259 sdrhoPmT =  NaN(1,ntn);
0260 ssdrhoPT =  NaN(1,ntn);
0261 srdrhoPmT =  NaN(1,ntn);
0262 %
0263 sJm =  NaN(ntn,nr);
0264 ssJ =  NaN(ntn,nr);
0265 srJm =  NaN(ntn,nr);
0266 %
0267 sPm =  NaN(ntn,nr);
0268 ssP =  NaN(ntn,nr);
0269 srPm =  NaN(ntn,nr);
0270 %
0271 for itn = 1:ntn,
0272     sI_totm(itn) = mean(sI_tot(1:itn));% time-integrated current average
0273     ssI_tot(itn) = std(sI_tot(1:itn));% time-integrated current STD
0274     srI_totm(itn) = ssI_tot(itn)/sI_totm(itn)/sqrt(itn);%relative statistical accuracy of current average
0275     %
0276     sP_rfm(itn) = mean(sP_rf(1:itn));% time-integrated power average
0277     ssP_rf(itn) = std(sP_rf(1:itn));% time-integrated power STD
0278     srP_rfm(itn) = ssP_rf(itn)/sP_rfm(itn)/sqrt(itn);%relative statistical accuracy of power average
0279     %
0280     srhoJm(itn) = mean(srhoJ(1:itn));% time-integrated current peak location average
0281     ssrhoJ(itn) = std(srhoJ(1:itn));% time-integrated current peak location STD
0282     srrhoJm(itn) = ssrhoJ(itn)/srhoJm(itn)/sqrt(itn);%relative statistical accuracy of current peak location average
0283     %
0284     srhoJmT(itn) = mean(srhoJT(1:itn));% time-integrated current peak location average
0285     ssrhoJT(itn) = std(srhoJT(1:itn));% time-integrated current peak location STD
0286     srrhoJmT(itn) = ssrhoJT(itn)/srhoJmT(itn)/sqrt(itn);%relative statistical accuracy of current peak location average
0287     %
0288     srhoPm(itn) = mean(srhoP(1:itn));% time-integrated power peak location average
0289     ssrhoP(itn) = std(srhoP(1:itn));% time-integrated power peak location STD
0290     srrhoPm(itn) = ssrhoP(itn)/srhoPm(itn)/sqrt(itn);%relative statistical accuracy of power peak location average
0291     %
0292     srhoPmT(itn) = mean(srhoPT(1:itn));% time-integrated power peak location average
0293     ssrhoPT(itn) = std(srhoPT(1:itn));% time-integrated power peak location STD
0294     srrhoPmT(itn) = ssrhoPT(itn)/srhoPmT(itn)/sqrt(itn);%relative statistical accuracy of power peak location average
0295     %
0296     sdrhoJm(itn) = mean(sdrhoJ(1:itn));% time-integrated current FWHM average
0297     ssdrhoJ(itn) = std(sdrhoJ(1:itn));% time-integrated current FWHM STD
0298     srdrhoJm(itn) = ssdrhoJ(itn)/sdrhoJm(itn)/sqrt(itn);%relative statistical accuracy of current FWHM average
0299     %
0300     sdrhoJmT(itn) = mean(sdrhoJT(1:itn));% time-integrated current FWHM average
0301     ssdrhoJT(itn) = std(sdrhoJT(1:itn));% time-integrated current FWHM STD
0302     srdrhoJmT(itn) = ssdrhoJT(itn)/sdrhoJmT(itn)/sqrt(itn);%relative statistical accuracy of current FWHM average
0303     %
0304     sdrhoPm(itn) = mean(sdrhoP(1:itn));% time-integrated power FWHM average
0305     ssdrhoP(itn) = std(sdrhoP(1:itn));% time-integrated power FWHM STD
0306     srdrhoPm(itn) = ssdrhoP(itn)/sdrhoPm(itn)/sqrt(itn);%relative statistical accuracy of power FWHM average
0307     %
0308     sdrhoPmT(itn) = mean(sdrhoPT(1:itn));% time-integrated power FWHM average
0309     ssdrhoPT(itn) = std(sdrhoPT(1:itn));% time-integrated power FWHM STD
0310     srdrhoPmT(itn) = ssdrhoPT(itn)/sdrhoPmT(itn)/sqrt(itn);%relative statistical accuracy of power FWHM average
0311     %
0312     sJm(itn,:) = mean(sJ(1:itn,:),1);% time-dependent average current profile
0313     ssJ(itn,:) = std(sJ(1:itn,:),0,1);% time-dependent average current STD
0314     srJm(itn,:) = ssJ(itn,:)./sJm(itn,:)/sqrt(itn);% relative statistical accuracy of average current profile
0315     %
0316     sPm(itn,:) = mean(sP(1:itn,:),1);% time-dependent average power profile
0317     ssP(itn,:) = std(sP(1:itn,:),0,1);% time-dependent average power STD
0318     srPm(itn,:) = ssP(itn,:)./sPm(itn,:)/sqrt(itn);% relative statistical accuracy of average power profile
0319     %
0320 end
0321 %
0322 sIm_tot = sJm*xdA;% time-dependent average current - should be the same as sI_totm
0323 sPm_rf = sPm*xdV;% time-dependent average power - should be the same as sP_rfm
0324 %
0325 smrhoJ = sJm*(xrho.*xdA)./sIm_tot;% time-dependent peak location of averaged current
0326 smrhoP = sPm*(xrho.*xdV)./sPm_rf;% time-dependent  peak location of averaged power
0327 %
0328 smrhoJT = sJm*(xrhoT.*xdA)./sIm_tot;% time-dependent peak location of averaged current
0329 smrhoPT = sPm*(xrhoT.*xdV)./sPm_rf;% time-dependent  peak location of averaged power
0330 %
0331 smxrhoJ = repmat(smrhoJ,[1,nr]);
0332 smxrhoP = repmat(smrhoP,[1,nr]);
0333 %
0334 smxrhoJT = repmat(smrhoJT,[1,nr]);
0335 smxrhoPT = repmat(smrhoPT,[1,nr]);
0336 %
0337 smdrhoJ = fwhm*sqrt((sJm.*((srho-smxrhoJ).^2)*xdA)./sIm_tot);%time-dependent FWHM of averaged current
0338 smdrhoP = fwhm*sqrt((sPm.*((srho-smxrhoP).^2)*xdV)./sPm_rf);%time-dependent FWHM of averaged power
0339 %
0340 smdrhoJT = fwhm*sqrt((sJm.*((srhoT-smxrhoJT).^2)*xdA)./sIm_tot);%time-dependent FWHM of averaged current
0341 smdrhoPT = fwhm*sqrt((sPm.*((srhoT-smxrhoPT).^2)*xdV)./sPm_rf);%time-dependent FWHM of averaged power
0342 %
0343 colors = {'k',[0,0.5,0]};
0344 %
0345 if fluct_display.mode >= 2,
0346     %
0347     % --------- current ---------
0348     %
0349     stnmask = 1:fix(ntn/10):ntn;
0350     %
0351     figure(110),clf
0352     %
0353     xlim = [0,ceil(max(stn))];
0354     xlab = 't/\tau_c';
0355     ylim = NaN;
0356     ylab = 'I (MA)';
0357     tit = '';
0358     leg = {'With fluct.','Average'};%{'Fluctuations'};%
0359     %
0360     leg_loc = leg;
0361     graph1D_jd(stn,sI_tot,0,0,xlab,ylab,tit,NaN,xlim,ylim,'-','+','r',2,20+14i,gca,0.9,0.7,0.7);
0362     %graph1D_jd(stn,sI_totm,0,0,'','','',NaN,xlim,ylim,'--','none','b',2,20+14i);
0363     %graph1D_jd(xlim,sI_totm(end)*[1,1],0,0,'','','',NaN,xlim,ylim,'--','none','b',2,20+14i);
0364     graph1D_jd(stn(stnmask),sI_totm(stnmask),0,0,'','','',leg_loc,xlim,ylim,'--','none','b',2,20+14i,gca,1,0.7,0.7,srI_totm(stnmask).*sI_totm(stnmask));
0365     if ~isempty(Zfluct_ref),
0366         leg_loc = [leg,leg_ref];
0367         for iref = 1:length(Zfluct_ref),
0368             graph1D_jd(xlim,Zfluct_ref(iref).J.sI_tot*[1,1],0,0,'','','',leg_loc,xlim,ylim,'-','none',colors{iref},2,20+14i);
0369         end
0370     else
0371         leg_loc = leg;
0372     end
0373     if fluct_display.stat,% to verify it is the same as sI_totm
0374         leg_loc = [leg_loc,'From av. profile']; 
0375         graph1D_jd(stn,sIm_tot,0,0,'','','',NaN,xlim,ylim,'--','none','m',0.5,20+14i);
0376     end
0377     %graph1D_jd(Ztime(Zmask(sitnPr)),I_tot(Zmask(sitnPr)),0,0,'','','',NaN,xlim,ylim,'none','o','m',2,20+14i);
0378     %graph1D_jd(Ztime(Zmask(sitnIr)),I_tot(Zmask(sitnIr)),0,0,'','','',NaN,xlim,ylim,'none','*','r',2,20+14i);
0379     %
0380     % --------- current radial location ---------
0381     %
0382     figure(120),clf
0383     %
0384     ylim = NaN;
0385     ylab = 'r/a';
0386     tit = 'current deposition';
0387     %
0388     graph1D_jd(stn,srhoJ,0,0,xlab,ylab,tit,NaN,xlim,ylim,'-','+','r',2,20+14i,gca,0.9,0.7,0.7);
0389     graph1D_jd(stn,srhoJm,0,0,'','','',NaN,xlim,ylim,'--','none','b',2,20+14i);
0390     %graph1D_jd(xlim,srhoJm(end)*[1,1],0,0,'','','',NaN,xlim,ylim,'--','none','b',2,20+14i);
0391     if ~isempty(Zfluct_ref),
0392         leg_loc = [leg,leg_ref];
0393         for iref = 1:length(Zfluct_ref),
0394             graph1D_jd(xlim,Zfluct_ref(iref).J.srho*[1,1],0,0,'','','',NaN,xlim,ylim,'-','none',colors{iref},2,20+14i);
0395         end        
0396     else
0397         leg_loc = leg;
0398     end
0399     if fluct_display.stat,% should be slightly different from srhoJm
0400         leg_loc = [leg_loc,'From av. profile']; 
0401         graph1D_jd(stn,smrhoJ,0,0,'','','',NaN,xlim,ylim,'--','none','m',0.5,20+14i);
0402     end
0403     graph1D_jd(stn(stnmask),srhoJm(stnmask),0,0,'','','',leg_loc,xlim,ylim,'none','none','b',2,20+14i,gca,1,0.7,0.7,srrhoJm(stnmask).*srhoJm(stnmask));
0404     %
0405     figure(121),clf
0406     %
0407     ylim = NaN;
0408     ylab = '\rho_{T}';
0409     tit = 'current deposition';
0410     %
0411     graph1D_jd(stn,srhoJT,0,0,xlab,ylab,tit,NaN,xlim,ylim,'-','+','r',2,20+14i,gca,0.9,0.7,0.7);
0412     graph1D_jd(stn,srhoJmT,0,0,'','','',NaN,xlim,ylim,'--','none','b',2,20+14i);
0413     leg_loc = leg;
0414     %graph1D_jd(stn(stnmask),srhoJmT(stnmask),0,0,'','','',leg_loc,xlim,ylim,'none','o','b',2,20+14i,gca,1,0.7,0.7,srrhoJmT(stnmask).*srhoJmT(stnmask));
0415     if ~isempty(Zfluct_ref),
0416         leg_loc = [leg,leg_ref];
0417         for iref = 1:length(Zfluct_ref),
0418             Zfluct_ref(iref).J.srhoT = ppval(Zfluct_ref(iref).J.srho,fitrhoT.p);
0419             graph1D_jd(xlim,Zfluct_ref(iref).J.srhoT*[1,1],0,0,'','','',leg_loc,xlim,ylim,'-','none',colors{iref},2,20+14i);
0420         end        
0421     end
0422     %
0423     % --------- current radial width ---------
0424     %
0425     figure(130),clf
0426     %
0427     ylim = NaN;
0428     ylab = 'FWHM (r/a)';
0429     %
0430     graph1D_jd(stn,sdrhoJ,0,0,xlab,ylab,tit,NaN,xlim,ylim,'-','+','r',2,20+14i,gca,0.9,0.7,0.7);
0431     graph1D_jd(stn,sdrhoJm,0,0,'','','',NaN,xlim,ylim,'--','none','b',2,20+14i);
0432     %graph1D_jd(xlim,sdrhoJm(end)*[1,1],0,0,'','','',NaN,xlim,ylim,'--','none','b',2,20+14i);
0433     if ~isempty(Zfluct_ref),
0434         leg_loc = [leg,leg_ref];
0435         for iref = 1:length(Zfluct_ref),
0436             graph1D_jd(xlim,Zfluct_ref(iref).J.sdrho*[1,1],0,0,'','','',NaN,xlim,ylim,'-','none',colors{iref},2,20+14i);
0437         end             
0438     else
0439         leg_loc = leg;
0440     end
0441     %
0442     if fluct_display.stat,
0443         leg_loc = [leg_loc,'\sigma_{\rho}','(\sigma_{\rho}^2 + FWHM^2)^{(1/2)}','From av. profile']; 
0444         graph1D_jd(stn,fwhm*ssrhoJ,0,0,'','','',NaN,xlim,ylim,'--','none','m',0.5,20+14i);
0445         graph1D_jd(stn,sqrt((fwhm*ssrhoJ).^2 + sdrhoJm.^2),0,0,'','','',NaN,xlim,ylim,'--','none','c',0.5,20+14i);
0446         graph1D_jd(stn,smdrhoJ,0,0,'','','',NaN,xlim,ylim,'-','none','y',2,20+14i);
0447     end
0448     graph1D_jd(stn(stnmask),sdrhoJm(stnmask),0,0,'','','',leg_loc,xlim,ylim,'none','none','b',2,20+14i,gca,1,0.7,0.7,srdrhoJm(stnmask).*sdrhoJm(stnmask));
0449     %
0450     figure(131),clf
0451     %
0452     ylim = NaN;
0453     ylab = 'FWHM (\rho_{T})';
0454     %
0455     graph1D_jd(stn,sdrhoJT,0,0,xlab,ylab,tit,NaN,xlim,ylim,'-','+','r',2,20+14i,gca,0.9,0.7,0.7);
0456     graph1D_jd(stn,sdrhoJmT,0,0,'','','',NaN,xlim,ylim,'--','none','b',2,20+14i);
0457     leg_loc = leg;
0458     %
0459     graph1D_jd(stn(stnmask),sdrhoJmT(stnmask),0,0,'','','',leg_loc,xlim,ylim,'none','none','b',2,20+14i,gca,1,0.7,0.7,srdrhoJmT(stnmask).*sdrhoJmT(stnmask));
0460     if ~isempty(Zfluct_ref),
0461         leg_loc = [leg,leg_ref];
0462         for iref = 1:length(Zfluct_ref),
0463             Zfluct_ref(iref).J.srhoJT = Zfluct_ref(iref).J.sJ*(Zfluct_ref(iref).rhoT'.*Zfluct_ref(iref).dA)./Zfluct_ref(iref).J.sI_tot;
0464             Zfluct_ref(iref).J.sdrhoT = fwhm*sqrt((Zfluct_ref(iref).J.sJ.*((Zfluct_ref(iref).rhoT-Zfluct_ref(iref).J.srhoJT).^2)*Zfluct_ref(iref).dA)./Zfluct_ref(iref).J.sI_tot);%time-dependent current FWHM rhoT grid
0465             graph1D_jd(xlim,Zfluct_ref(iref).J.sdrhoT*[1,1],0,0,'','','',NaN,xlim,ylim,'-','none',colors{iref},2,20+14i);
0466         end             
0467     else
0468         leg_loc = leg;
0469     end
0470     %
0471     % --------- power ---------
0472     %
0473     figure(140),clf
0474     %
0475     ylim = NaN;
0476     ylab = 'P_{RF} (MW)';
0477     tit = '';
0478     %
0479     graph1D_jd(stn,sP_rf,0,0,xlab,ylab,tit,NaN,xlim,ylim,'-','+','r',2,20+14i,gca,0.9,0.7,0.7);
0480     graph1D_jd(stn,sP_rfm,0,0,'','','',NaN,xlim,ylim,'--','none','b',2,20+14i);
0481     %graph1D_jd(xlim,sP_rfm(end)*[1,1],0,0,'','','',NaN,xlim,ylim,'--','none','b',2,20+14i);
0482     if ~isempty(Zfluct_ref),
0483         leg_loc = [leg,leg_ref];
0484         for iref = 1:length(Zfluct_ref),
0485             graph1D_jd(xlim,Zfluct_ref(iref).P.sP_rf*[1,1],0,0,'','','',NaN,xlim,ylim,'-','none',colors{iref},2,20+14i);
0486         end             
0487     else
0488         leg_loc = leg;
0489     end
0490     if fluct_display.stat,% to verify it is the same as sP_rfm
0491         leg_loc = [leg_loc,'From av. profile']; 
0492         graph1D_jd(stn,sPm_rf,0,0,'','','',NaN,xlim,ylim,'--','none','m',0.5,20+14i);
0493     end
0494     graph1D_jd(stn(stnmask),sP_rfm(stnmask),0,0,'','','',leg_loc,xlim,ylim,'none','none','b',2,20+14i,gca,1,0.7,0.7,srP_rfm(stnmask).*sP_rfm(stnmask));
0495     %graph1D_jd(Ztime(Zmask(sitnPr)),P_rf(Zmask(sitnPr)),0,0,'','','',NaN,xlim,ylim,'none','o','r',2,20+14i);
0496     %graph1D_jd(Ztime(Zmask(sitnIr)),P_rf(Zmask(sitnIr)),0,0,'','','',NaN,xlim,ylim,'none','*','m',2,20+14i);
0497     %
0498     % --------- power radial location ---------
0499     %
0500     figure(150),clf
0501     %
0502     ylim = NaN;
0503     ylab = 'r/a';
0504     tit = 'power deposition';
0505     %
0506     graph1D_jd(stn,srhoP,0,0,xlab,ylab,tit,NaN,xlim,ylim,'-','+','r',2,20+14i,gca,0.9,0.7,0.7);
0507     graph1D_jd(stn,srhoPm,0,0,'','','',NaN,xlim,ylim,'--','none','b',2,20+14i);
0508     %graph1D_jd(xlim,srhoPm(end)*[1,1],0,0,'','','',NaN,xlim,ylim,'--','none','b',2,20+14i);
0509     if ~isempty(Zfluct_ref),
0510         leg_loc = [leg,leg_ref];
0511         for iref = 1:length(Zfluct_ref),
0512             graph1D_jd(xlim,Zfluct_ref(iref).P.srho*[1,1],0,0,'','','',NaN,xlim,ylim,'-','none',colors{iref},2,20+14i);
0513         end             
0514     else
0515         leg_loc = leg;
0516     end
0517     %
0518     if fluct_display.stat,% should be slightly different from srhoPm
0519         leg_loc = [leg_loc,'From av. profile']; 
0520         graph1D_jd(stn,smrhoP,0,0,'','','',NaN,xlim,ylim,'--','none','m',0.5,20+14i);
0521     end
0522     %
0523     graph1D_jd(stn(stnmask),srhoPm(stnmask),0,0,'','','',leg_loc,xlim,ylim,'none','none','b',2,20+14i,gca,1,0.7,0.7,srrhoPm(stnmask).*srhoPm(stnmask));
0524     %
0525     figure(151),clf
0526     %
0527     ylim = NaN;
0528     ylab = '\rho_{T}';
0529     tit = 'power deposition';
0530     %
0531     graph1D_jd(stn,srhoPT,0,0,xlab,ylab,tit,NaN,xlim,ylim,'-','+','r',2,20+14i,gca,0.9,0.7,0.7);
0532     graph1D_jd(stn,srhoPmT,0,0,'','','',NaN,xlim,ylim,'--','none','b',2,20+14i);
0533     leg_loc = leg;
0534     %
0535     if ~isempty(Zfluct_ref),
0536         leg_loc = [leg,leg_ref];
0537         for iref = 1:length(Zfluct_ref),
0538             Zfluct_ref(iref).P.srhoT = ppval(Zfluct_ref(iref).P.srho,fitrhoT.p);
0539             graph1D_jd(xlim,Zfluct_ref(iref).P.srhoT*[1,1],0,0,'','','',leg_loc,xlim,ylim,'-','none',colors{iref},2,20+14i);
0540         end             
0541     else
0542         leg_loc = leg;
0543     end
0544     graph1D_jd(stn(stnmask),srhoPmT(stnmask),0,0,'','','',leg_loc,xlim,ylim,'none','none','b',2,20+14i,gca,1,0.7,0.7,srrhoPmT(stnmask).*srhoPmT(stnmask));
0545 
0546     %
0547     % --------- power radial width ---------
0548     %
0549     figure(160),clf
0550     %
0551     ylim = NaN;
0552     ylab = 'FWHM (r/a)';
0553     %
0554     graph1D_jd(stn,sdrhoP,0,0,xlab,ylab,tit,NaN,xlim,ylim,'-','+','r',2,20+14i,gca,0.9,0.7,0.7);
0555     graph1D_jd(stn,sdrhoPm,0,0,'','','',NaN,xlim,ylim,'--','none','b',2,20+14i);
0556     %graph1D_jd(xlim,sdrhoPm(end)*[1,1],0,0,'','','',NaN,xlim,ylim,'--','none','b',2,20+14i);
0557     if ~isempty(Zfluct_ref),
0558         leg_loc = [leg,leg_ref];
0559         for iref = 1:length(Zfluct_ref),
0560             graph1D_jd(xlim,Zfluct_ref(iref).P.sdrho*[1,1],0,0,'','','',NaN,xlim,ylim,'-','none',colors{iref},2,20+14i);
0561         end             
0562     else
0563         leg_loc = leg;
0564     end
0565     if fluct_display.stat,
0566         leg_loc = [leg_loc,'\sigma_{\rho}','(\sigma_{\rho}^2 + FWHM^2)^{(1/2)}','From av. profile']; 
0567         graph1D_jd(stn,fwhm*ssrhoP,0,0,'','','',NaN,xlim,ylim,'--','none','m',0.5,20+14i);
0568         graph1D_jd(stn,sqrt((fwhm*ssrhoP).^2 + sdrhoPm.^2),0,0,'','','',NaN,xlim,ylim,'--','none','c',0.5,20+14i);
0569         graph1D_jd(stn,smdrhoP,0,0,'','','',NaN,xlim,ylim,'-','none','y',2,20+14i);
0570     end
0571     graph1D_jd(stn(stnmask),sdrhoPm(stnmask),0,0,'','','',leg_loc,xlim,ylim,'none','none','b',2,20+14i,gca,1,0.7,0.7,srdrhoPm(stnmask).*sdrhoPm(stnmask));
0572     %
0573     figure(161),clf
0574     %
0575     ylim = NaN;
0576     ylab = 'FWHM (\rho_{T})';
0577     %
0578     graph1D_jd(stn,sdrhoPT,0,0,xlab,ylab,tit,NaN,xlim,ylim,'-','+','r',2,20+14i,gca,0.9,0.7,0.7);
0579     graph1D_jd(stn,sdrhoPmT,0,0,'','','',NaN,xlim,ylim,'--','none','b',2,20+14i);
0580     leg_loc = leg;
0581     %
0582     if ~isempty(Zfluct_ref),
0583         leg_loc = [leg,leg_ref];
0584         for iref = 1:length(Zfluct_ref),
0585             srhoP = sP*(xrho.*xdV)./sP_rf;% time-dependent power peak location
0586             Zfluct_ref(iref).P.srhoPT = Zfluct_ref(iref).P.sP*(Zfluct_ref(iref).rhoT'.*Zfluct_ref(iref).dV)./Zfluct_ref(iref).P.sP_rf;
0587             Zfluct_ref(iref).P.sdrhoT = fwhm*sqrt((Zfluct_ref(iref).P.sP.*((Zfluct_ref(iref).rhoT-Zfluct_ref(iref).P.srhoPT).^2)*Zfluct_ref(iref).dV)./Zfluct_ref(iref).P.sP_rf);%time-dependent power FWHM rhoT grid
0588             graph1D_jd(xlim,Zfluct_ref(iref).P.sdrhoT*[1,1],0,0,'','','',leg_loc,xlim,ylim,'-','none',colors{iref},2,20+14i);
0589         end             
0590     else
0591         leg_loc = leg;
0592     end
0593     %
0594     % --------- current profile ---------
0595     %
0596     figure(170),clf
0597     %
0598     xlim = [0,1];
0599     xlab = 'r/a';
0600     ylim = NaN;
0601     ylab = 'dJ/d\rho (MA)';
0602     tit = '';
0603     leg = {'With fluct.'};
0604     %
0605     sxdAdrho = repmat((xdA./xdrho).',[ntn,1]);
0606     sxdIdrhom = sJm.*sxdAdrho;
0607     %
0608     smxdrhoJ = repmat(smdrhoJ,[1,nr]);
0609     sxIm_tot = repmat(sIm_tot,[1,nr]);
0610     %
0611     sxdIdrhom_fit = sxIm_tot.*exp(-(srho - smxrhoJ).^2./(smxdrhoJ*sqrt(2)/fwhm).^2)./(sqrt(pi)*smxdrhoJ*sqrt(2)/fwhm);
0612     %
0613     %graph1D_jd(xrho,sJm(end,:),0,0,xlab,ylab,tit,NaN,xlim,ylim,'-','+','r',2,20+14i,gca,0.9,0.7,0.7,ssJ(end,:));
0614     graph1D_jd(xrho,sxdIdrhom(end,:),0,0,xlab,ylab,tit,NaN,xlim,ylim,'-','+','r',2,20+14i,gca,0.9,0.7,0.7,sxdIdrhom(end,:).*srJm(end,:));
0615     if ~isempty(Zfluct_ref),
0616         leg_loc = [leg,leg_ref];
0617         for iref = 1:length(Zfluct_ref),
0618             xdIdrhom_ref = Zfluct_ref(iref).J.sJm.*(Zfluct_ref(iref).dA./Zfluct_ref(iref).drho).';
0619             graph1D_jd(Zfluct_ref(iref).rho,xdIdrhom_ref,0,0,'','','',NaN,xlim,ylim,'-','none',colors{iref},2,20+14i);
0620         end             
0621     else
0622         leg_loc = leg;
0623     end
0624     if fluct_display.stat,
0625         leg_loc = [leg_loc,'Fit']; 
0626         graph1D_jd(xrho,sxdIdrhom_fit(end,:),0,0,'','','',leg_loc,xlim,ylim,'--','none','b',2,20+14i);
0627     else
0628         graph1D_jd(NaN,NaN,0,0,'','','',leg_loc,xlim,ylim,'--','none','b',2,20+14i);
0629     end
0630     %
0631     figure(171),clf
0632     %
0633     xlab = '\rho_{T}';
0634     ylab = 'dJ/d\rho_{T} (MA)';
0635     %
0636     sxdAdrhoT = repmat((xdA./xdrhoT).',[ntn,1]);
0637     sxdIdrhomT = sJm.*sxdAdrhoT;
0638     %
0639     smxdrhoJT = repmat(smdrhoJT,[1,nr]);
0640     sxIm_tot = repmat(sIm_tot,[1,nr]);
0641     %
0642     sxdIdrhom_fitT = sxIm_tot.*exp(-(srhoT - smxrhoJT).^2./(smxdrhoJT*sqrt(2)/fwhm).^2)./(sqrt(pi)*smxdrhoJT*sqrt(2)/fwhm);
0643     %
0644     graph1D_jd(xrhoT,sxdIdrhomT(end,:),0,0,xlab,ylab,tit,NaN,xlim,ylim,'-','+','r',2,20+14i,gca,0.9,0.7,0.7,sxdIdrhomT(end,:).*srJm(end,:));
0645     leg_loc = leg;  
0646     %
0647     if ~isempty(Zfluct_ref),
0648         leg_loc = [leg,leg_ref];
0649         for iref = 1:length(Zfluct_ref),
0650             Zfluct_ref(iref).rhoT_S = [0,interp1(xrhoT',xrhoT_S(2:end),Zfluct_ref(iref).rhoT,'spline')];
0651             Zfluct_ref(iref).drhoT = diff(Zfluct_ref(iref).rhoT_S.');       
0652             xdIdrhom_refT = Zfluct_ref(iref).J.sJm.*(Zfluct_ref(iref).dA./Zfluct_ref(iref).drhoT).';
0653             graph1D_jd(Zfluct_ref(iref).rhoT,xdIdrhom_refT,0,0,'','','',leg_loc,xlim,ylim,'-','none',colors{iref},2,20+14i);
0654         end             
0655     else
0656         leg_loc = leg;
0657     end
0658     graph1D_jd(NaN,NaN,0,0,'','','',leg_loc,xlim,ylim,'--','none','b',2,20+14i);
0659 
0660     %
0661     % --------- power profile ---------
0662     %
0663     figure(180),clf
0664     %
0665     xlim = [0,1];
0666     xlab = 'r/a';
0667     ylim = NaN;
0668     ylab = 'dP_{RF}/d\rho (MW)';
0669     tit = '';
0670     leg = {'With fluct.'};
0671     %
0672     sxdVdrho = repmat((xdV./xdrho).',[ntn,1]);
0673     sxdPdrhom = sPm.*sxdVdrho;
0674     %
0675     smxdrhoP = repmat(smdrhoP,[1,nr]);
0676     sxPm_rf = repmat(sPm_rf,[1,nr]);
0677     %
0678     sxdPdrhom_fit = sxPm_rf.*exp(-(srho - smxrhoP).^2./(smxdrhoP*sqrt(2)/fwhm).^2)./(sqrt(pi)*smxdrhoP*sqrt(2)/fwhm);
0679     %
0680     %graph1D_jd(xrho,sJm(end,:),0,0,xlab,ylab,tit,NaN,xlim,ylim,'-','+','r',2,20+14i,gca,0.9,0.7,0.7,ssJ(end,:));
0681     graph1D_jd(xrho,sxdPdrhom(end,:),0,0,xlab,ylab,tit,NaN,xlim,ylim,'-','+','r',2,20+14i,gca,0.9,0.7,0.7,sxdPdrhom(end,:).*srPm(end,:));
0682     if ~isempty(Zfluct_ref),
0683         leg_loc = [leg,leg_ref];
0684         for iref = 1:length(Zfluct_ref),
0685             xdPdrhom_ref = Zfluct_ref(iref).P.sPm.*(Zfluct_ref(iref).dV./Zfluct_ref(iref).drho).';
0686             graph1D_jd(Zfluct_ref(iref).rho,xdPdrhom_ref,0,0,'','','',NaN,xlim,ylim,'-','none',colors{iref},2,20+14i);
0687         end             
0688     else
0689         leg_loc = leg;
0690     end
0691     if fluct_display.stat,
0692         leg_loc = [leg_loc,'Fit']; 
0693         graph1D_jd(xrho,sxdPdrhom_fit(end,:),0,0,'','','',leg_loc,xlim,ylim,'--','none','b',2,20+14i);
0694     else
0695         graph1D_jd(NaN,NaN,0,0,'','','',leg_loc,xlim,ylim,'--','none','b',2,20+14i);
0696     end
0697     %
0698     figure(181),clf
0699     %
0700     xlab = '\rho_{T}';
0701     ylab = 'dP_{RF}/d\rho_{T} (MW)';
0702     %
0703     sxdVdrhoT = repmat((xdV./xdrhoT).',[ntn,1]);
0704     sxdPdrhomT = sPm.*sxdVdrhoT;
0705     %
0706     smxdrhoPT = repmat(smdrhoPT,[1,nr]);
0707     sxPm_rf = repmat(sPm_rf,[1,nr]);
0708     %
0709     sxdPdrhom_fitT = sxPm_rf.*exp(-(srhoT - smxrhoPT).^2./(smxdrhoPT*sqrt(2)/fwhm).^2)./(sqrt(pi)*smxdrhoPT*sqrt(2)/fwhm);
0710     %
0711     graph1D_jd(xrhoT,sxdPdrhomT(end,:),0,0,xlab,ylab,tit,NaN,xlim,ylim,'-','+','r',2,20+14i,gca,0.9,0.7,0.7,sxdPdrhomT(end,:).*srPm(end,:));
0712     if ~isempty(Zfluct_ref),
0713         leg_loc = [leg,leg_ref];
0714         for iref = 1:length(Zfluct_ref),
0715             xdPdrhom_refT = Zfluct_ref(iref).P.sPm.*(Zfluct_ref(iref).dV./Zfluct_ref(iref).drhoT).';
0716             graph1D_jd(Zfluct_ref(iref).rhoT,xdPdrhom_refT,0,0,'','','',leg_loc,xlim,ylim,'-','none',colors{iref},2,20+14i);
0717         end             
0718     else
0719         leg_loc = leg;
0720     end
0721     %
0722     %
0723     % --------- power density profile ---------
0724     %
0725     figure(190),clf
0726     %
0727     xlim = [0,1];
0728     xlab = 'r/a';
0729     ylim = NaN;
0730     ylab = 'P_{RF} (MW/m^3)';
0731     tit = '';
0732     %
0733     if isfield(opt,'hxrinv');
0734         %
0735         hxrinv = load(opt.hxrinv);
0736         hxrinv.rho = hxrinv.rayon/(1000*hxrinv.ap);
0737         hxrinv.Xinterp = interp1(hxrinv.rho,hxrinv.X,xrho).';
0738         hxrinv.P = hxrinv.Xinterp*(sPm(end,:)*xdrho)/(hxrinv.Xinterp*xdrho);
0739         hxrinv.P(find(hxrinv.P < 0,1,'first'):end) = 0;
0740         %
0741         graph1D_jd(xrho,hxrinv.P,0,0,'','','',NaN,xlim,ylim,'-','none','b',2,20+14i);
0742         %
0743         leg_loc = {'Exp. HXR'};
0744         %
0745     else
0746         leg_loc = {};
0747     end
0748     %
0749     if ~isempty(Zfluct_ref),
0750         leg_loc = [leg_loc,leg_ref];
0751         for iref = 1:length(Zfluct_ref),
0752             graph1D_jd(Zfluct_ref(iref).rho,Zfluct_ref(iref).P.sPm,0,0,'','','',NaN,xlim,ylim,'-','none',colors{iref},2,20+14i);
0753         end             
0754     end
0755     %
0756     leg = [leg_loc,'With fluct.'];
0757     graph1D_jd(xrho,sPm(end,:),0,0,xlab,ylab,tit,leg,xlim,ylim,'-','+','r',2,20+14i,gca,0.9,0.7,0.7,sxdPdrhom(end,:).*srPm(end,:));
0758     %
0759     figure(191),clf
0760     %
0761     xlab = '\rho_{T}';
0762     ylab = 'P_{RF} (MW/m^3)';
0763     tit = '';
0764     %
0765     leg_loc = {};
0766     %
0767     if ~isempty(Zfluct_ref),
0768         leg_loc = [leg_loc,leg_ref];
0769         for iref = 1:length(Zfluct_ref),
0770             graph1D_jd(Zfluct_ref(iref).rhoT,Zfluct_ref(iref).P.sPm,0,0,'','','',NaN,xlim,ylim,'-','none',colors{iref},2,20+14i);
0771         end             
0772     end
0773     %
0774     leg = [leg_loc,'With fluct.'];
0775     graph1D_jd(xrhoT,sPm(end,:),0,0,xlab,ylab,tit,leg,xlim,ylim,'-','+','r',2,20+14i,gca,0.9,0.7,0.7,sxdPdrhom(end,:).*srPm(end,:));
0776     %
0777 end
0778 %
0779 Zfluct_out.stn = stn;
0780 Zfluct_out.rho = xrho;
0781 Zfluct_out.drho = xdrho;
0782 if isfield(equilDKE,'xrhoT')
0783     Zfluct_out.rhoT = equilDKE.xrhoT;
0784 end
0785 Zfluct_out.dA = xdA;
0786 Zfluct_out.dV = xdV;
0787 Zfluct_out.calctime = calctime;
0788 Zfluct_out.J.sJ = sJ;
0789 Zfluct_out.J.sJm = sJm;
0790 Zfluct_out.J.ssJ = ssJ;
0791 Zfluct_out.J.srJm = srJm;
0792 Zfluct_out.J.sI_tot = sI_tot;
0793 Zfluct_out.J.sI_totm = sI_totm;
0794 Zfluct_out.J.ssI_tot = ssI_tot;
0795 Zfluct_out.J.srI_totm = srI_totm;
0796 Zfluct_out.J.srho = srhoJ;
0797 Zfluct_out.J.sdrho = sdrhoJ;
0798 Zfluct_out.J.srhom = srhoJm;
0799 Zfluct_out.J.ssrho = ssrhoJ;
0800 Zfluct_out.J.srrhom = srrhoJm;
0801 Zfluct_out.J.sdrhom = sdrhoJm;
0802 Zfluct_out.J.ssdrho = ssdrhoJ;
0803 Zfluct_out.J.srdrhom = srdrhoJm;
0804 Zfluct_out.J.smrho = smrhoJ;
0805 Zfluct_out.J.smdrho = smdrhoJ;
0806 Zfluct_out.P.sP = sP;
0807 Zfluct_out.P.sPm = sPm;
0808 Zfluct_out.P.ssP = ssP;
0809 Zfluct_out.P.srPm = srPm;
0810 Zfluct_out.P.sP_rf = sP_rf;
0811 Zfluct_out.P.sP_rfm = sP_rfm;
0812 Zfluct_out.P.ssP_rf = ssP_rf;
0813 Zfluct_out.P.srP_rfm = srP_rfm;
0814 Zfluct_out.P.srho = srhoP;
0815 Zfluct_out.P.sdrho = sdrhoP;
0816 Zfluct_out.P.srhom = srhoPm;
0817 Zfluct_out.P.ssrho = ssrhoP;
0818 Zfluct_out.P.srrhom = srrhoPm;
0819 Zfluct_out.P.sdrhom = sdrhoPm;
0820 Zfluct_out.P.ssdrho = ssdrhoP;
0821 Zfluct_out.P.srdrhom = srdrhoPm;
0822 Zfluct_out.P.smrho = smrhoP;
0823 Zfluct_out.P.smdrho = smdrhoP;
0824 %
0825 
0826