wave_disp_outputs_mt

0001 function wave_disp_outputs_mt(wavess,shotimes,style,axs)
0002 %
0003 nt = length(shotimes);
0004 %
0005 nr_dep = 40;
0006 %
0007 opt_peak = 0;
0008 %
0009 rho_dep_S = linspace(0,1,nr_dep+1);
0010 rho_dep = (rho_dep_S(1:end-1) + rho_dep_S(2:end))/2;
0011 drho_dep = diff(rho_dep_S);
0012 %
0013 alpha_L = NaN(nt,0);
0014 beta_L = NaN(nt,0);
0015 P0 = NaN(nt,0);
0016 Pabs = NaN(nt,0);
0017 rhoabs = NaN(nt,0);
0018 %
0019 launchids = {};
0020 %
0021 for it = 1:nt,
0022     waves = wavess{it};
0023     if isempty(waves),% times for which no RF wave is calculated
0024         %
0025         continue
0026         %
0027     else% times for which RF wave are calculated
0028         %
0029         for ilaun = 1:length(waves),
0030             %
0031             wave = waves{ilaun};
0032             launch = wave.rayinit.launch;
0033             %
0034             launchid = launch.id;
0035             %launchid = launchid(1:min([end,strfind(launchid,'_')-1]));
0036             %
0037             iw = find(strcmp(launchids,launchid));
0038             %
0039             if isempty(iw),% add new launcher to list
0040                 %
0041                 launchids = [launchids,launchid];
0042                 %
0043                 alpha_L(:,end+1) = NaN(nt,1);
0044                 beta_L(:,end+1) = NaN(nt,1);
0045                 P0(:,end+1) = NaN(nt,1);
0046                 Pabs(:,end+1) = NaN(nt,1);
0047                 rhoabs(:,end+1) = NaN(nt,1);
0048                 %
0049                 iw = length(launchids);
0050                 %
0051             elseif P0(it,iw) > 0
0052                 error('problem identifying launchers')
0053             end
0054             %
0055             alpha_L(it,iw) = launch.yalpha_L;
0056             beta_L(it,iw) = launch.ybeta_L;
0057             %
0058             dP_drho_lin = zeros(1,nr_dep);
0059             P0(it,iw) = 0;
0060             Pleft = 0;
0061             %
0062             if wave.rayinit.yP0_2piRp > 0
0063                 Rp2pi = launch.yP_L/wave.rayinit.yP0_2piRp;
0064                 %
0065                 for iy = 1:length(wave.rays),
0066                     P0(it,iw) = P0(it,iw) + Rp2pi*wave.rays{iy}.P0_2piRp;
0067                     srho_P = sqrt(wave.rays{iy}.spsin);
0068                     sP_norm = wave.rays{iy}.sP_2piRp_lin/wave.rays{iy}.sP_2piRp_lin(1);
0069                     dP_drho_lin = dP_drho_lin + Rp2pi*wave.rays{iy}.P0_2piRp*radialdampingprofile_jd(srho_P,sP_norm,rho_dep_S);
0070                     Pleft = Pleft + Rp2pi*wave.rays{iy}.sP_2piRp_lin(end);
0071                 end
0072                 %
0073                 Pabs(it,iw) = sum(dP_drho_lin.*drho_dep);
0074                 if abs(P0(it,iw) - Pabs(it,iw) - Pleft)/(P0(it,iw) + eps) > 1e-6,%to check consistency
0075                     disp('WARNING : consistency pb in power deposition')
0076                     keyboard
0077                 end
0078                 %
0079                 if all(isnan(dP_drho_lin)),
0080                     rhoabs(it,iw) = NaN;
0081                 elseif opt_peak == 0,%max position
0082                     rhoabs(it,iw) = rho_dep(find(dP_drho_lin == max(dP_drho_lin),1,'first'));
0083                 else
0084                     rhoabs(it,iw) = sum(dP_drho_lin.*rho_dep.*drho_dep)/Pabs(it,iw);
0085                 end
0086             else
0087                 rhoabs(it,iw) = NaN;
0088             end
0089         end
0090     end
0091 end
0092 %
0093 nw = length(launchids);
0094 %
0095 phi_L = pi - alpha_L;
0096 phi_L(phi_L > pi) = phi_L(phi_L > pi) - 2*pi;
0097 theta_L = pi/2 - beta_L;
0098 %
0099 fabs = Pabs./P0;
0100 %
0101 colors = {'r','b',[0,0.5,0],'k'};
0102 markers = {'s','d','+','o'};
0103 %
0104 xlab = 't (s)';
0105 xlim = 0;
0106 %
0107 if nargin < 4,
0108     red = 0.9;
0109     lspace = 0.7;
0110     lspace2 = 0.5;
0111     bspace = 0.7;
0112     bspace2 = 0.5;
0113 else
0114     red = 1;
0115     lspace = NaN;
0116     lspace2 = NaN;
0117     bspace = NaN;
0118     bspace2 = NaN;
0119 end
0120 %
0121 % poloidal angle
0122 %
0123 if nargin < 4,
0124     figure(1),clf,set(1,'Name','Poloidal angle')
0125     ax = gca;
0126 else
0127     ax = axs(1);
0128 end
0129 %
0130 ylab = 'pol. angle (deg.)';
0131 ylim = NaN;%[-90,90];
0132 %
0133 for iw = 1:nw-1,
0134     [ax] = graph1D_jd(shotimes,theta_L(:,iw)*180/pi,0,0,'','','',NaN,NaN,NaN,'-',markers{iw},colors{iw},2,style,ax);
0135 end
0136 [ax] = graph1D_jd(shotimes,theta_L(:,nw)*180/pi,0,0,xlab,ylab,'',launchids,xlim,ylim,'-',markers{nw},colors{nw},2,style,ax);
0137 %
0138 % toroidal angle
0139 %
0140 if nargin < 4,
0141     figure(2),clf,set(2,'Name','Toroidal angle')
0142     ax = gca;
0143 else
0144     ax = axs(2);
0145 end
0146 %
0147 ylab = 'tor. angle (deg.)';
0148 ylim = NaN;%[-90,90];
0149 %
0150 for iw = 1:nw-1,
0151     [ax] = graph1D_jd(shotimes,phi_L(:,iw)*180/pi,0,0,'','','',NaN,NaN,NaN,'-',markers{iw},colors{iw},2,style,ax);
0152 end
0153 [ax] = graph1D_jd(shotimes,phi_L(:,nw)*180/pi,0,0,xlab,ylab,'',launchids,xlim,ylim,'-',markers{nw},colors{nw},2,style,ax);
0154 %
0155 % launched power
0156 %
0157 if nargin < 4,
0158     figure(3),clf,set(3,'Name','launched power')
0159     ax = gca;
0160 else
0161     ax = axs(3);
0162 end
0163 %
0164 ylab = 'RF power (kW)';
0165 ylim = 0;
0166 %
0167 for iw = 1:nw-1,
0168     [ax] = graph1D_jd(shotimes,P0(:,iw)/1e3,0,0,'','','',NaN,NaN,NaN,'-',markers{iw},colors{iw},2,style,ax);
0169 end
0170 [ax] = graph1D_jd(shotimes,P0(:,nw)/1e3,0,0,xlab,ylab,'',launchids,xlim,ylim,'-',markers{nw},colors{nw},2,style,ax);
0171 %
0172 % power absorption
0173 %
0174 if nargin < 4,
0175     figure(4),clf,set(4,'Name','power absorption')
0176     ax = gca;
0177 else
0178     ax = axs(4);
0179 end
0180 %
0181 ylab = 'Power absorption (%)';
0182 ylim = [-10,110];
0183 %
0184 for iw = 1:nw-1,
0185     [ax] = graph1D_jd(shotimes,fabs(:,iw)*100,0,0,'','','',NaN,NaN,NaN,'-',markers{iw},colors{iw},2,style,ax);
0186 end
0187 [ax] = graph1D_jd(shotimes,fabs(:,nw)*100,0,0,xlab,ylab,'',launchids,xlim,ylim,'-',markers{nw},colors{nw},2,style,ax);
0188 %
0189 % power deposition
0190 %
0191 if nargin < 4,
0192     figure(5),clf,set(5,'Name','power deposition')
0193     ax = gca;
0194 else
0195     ax = axs(5);
0196 end
0197 %
0198 ylab = '\rho_{abs} (pol. flux)';
0199 ylim = [0,1];
0200 %
0201 for iw = 1:nw-1,
0202     [ax] = graph1D_jd(shotimes,rhoabs(:,iw),0,0,'','','',NaN,NaN,NaN,'-',markers{iw},colors{iw},2,style,ax);
0203 end
0204 [ax] = graph1D_jd(shotimes,rhoabs(:,nw),0,0,xlab,ylab,'',launchids,xlim,ylim,'-',markers{nw},colors{nw},2,style,ax);
0205 %