alphaphi_luke_jd

 This function uses the LUKE power deposition output to calculate the power deposited along the rays

0001 function waves = alphaphi_luke_jd(waves,dke_out,ZP0,equilDKE,mksa)
0002 %
0003 % This function uses the LUKE power deposition output to calculate the power deposited along the rays
0004 %
0005 nw = length(waves);
0006 if isempty(waves),
0007     return
0008 end
0009 %
0010 if ~isfield(dke_out,'zS'),
0011     waves = {};;
0012     return;
0013 end
0014 %
0015 zS = dke_out.zS;
0016 residue_rf = dke_out.residue_rf;
0017 xyprop_dke = dke_out.xyprop_dke;
0018 yb = dke_out.yb;
0019 nit = size(dke_out.zP_0_2piRp_mod,2);
0020 ny = size(xyprop_dke,2);
0021 nb = length(yb);
0022 
0023 %
0024 xydV_2piRp_dke = equilDKE.xdV_2piRp_dke.'*ones(1,ny);
0025 %
0026 for itn = 1:nit,
0027     %
0028     xyPabs_rf_dke_fsav_mksa = ZP0(itn).xy_rf_fsav*mksa.P_ref*1e6;% absorbed power density in W/m3 (RLA)
0029     xyPrf_2piRp = xyPabs_rf_dke_fsav_mksa.*xydV_2piRp_dke;%power absorbed from each ray within each FS (RLA)
0030     xyPrf_2piRp(xyPrf_2piRp < 0) = 0;% RLA
0031     %
0032     if iscell(residue_rf),
0033         nit_rf = length(residue_rf{itn});
0034     else
0035         nit_rf  = length(residue_rf);% old simulation
0036     end
0037     %
0038     if size(dke_out.zP_0_2piRp_mod,1) >= nit_rf,% all RF states kept in zP_0_2piRp_mod (RLA)
0039         zP_0_2piRp_mod = dke_out.zP_0_2piRp_mod{nit_rf,itn};% (RLA)
0040         if isfield(dke_out,'zP_0_2piRp_mod_coll'),
0041             zP_0_2piRp_mod_coll = dke_out.zP_0_2piRp_mod_coll{nit_rf,itn};% (RLA + NRCA)
0042         end
0043     else
0044         zP_0_2piRp_mod = dke_out.zP_0_2piRp_mod{1,itn};% old case where only last  RF state kept in zP_0_2piRp_mod (RLA)
0045         if isfield(dke_out,'zP_0_2piRp_mod_coll'),
0046             zP_0_2piRp_mod_coll = dke_out.zP_0_2piRp_mod_coll{1,itn};% old case where only last  RF state kept in zP_0_2piRp_mod (RLA + NRCA)
0047         end
0048     end
0049     %
0050     ibb = 0;
0051     %
0052     for iw = 1:nw,
0053         for ib = 1:length(waves{iw}.rays),
0054             ibb = ibb + 1;
0055             ray = waves{iw}.rays{ib};
0056             ns = length(ray.ss);
0057             %
0058             if ibb < nb,
0059                 iy_end = yb(ibb + 1) - 1;
0060             else
0061                 iy_end = ny;
0062             end
0063             %
0064             if any(xyprop_dke(:,iy_end) == 1)%outward propagating ray fragment (RLA)
0065                 Prf_end_2piRp = xyPrf_2piRp(find(xyprop_dke(:,iy_end) == 1,1,'last'),iy_end); 
0066             else%inward propagating ray fragment (RLA)
0067                 Prf_end_2piRp = xyPrf_2piRp(find(xyprop_dke(:,iy_end) == -1,1,'first'),iy_end);
0068             end
0069             %
0070             zP_end_2piRp = max([0,zP_0_2piRp_mod{ibb}(end) - Prf_end_2piRp/2]);%to account for all QL absorbed power (RLA)
0071             %
0072             zS_fit = [ray.ss(1),zS{ibb},ray.ss(end)];
0073             zP_0_2piRp_fit = [ray.P0_2piRp,zP_0_2piRp_mod{ibb},zP_end_2piRp];% RLA
0074             %
0075             % power along the ray
0076             %
0077             sP_2piRp = interp1(zS_fit,zP_0_2piRp_fit,waves{iw}.rays{ib}.ss,'pchip');% RLA
0078             %
0079             % optical depth
0080             %
0081             mask_ninf = sP_2piRp > 0;
0082             stau = Inf(1,ns);
0083             stau(mask_ninf) = log(ray.P0_2piRp./sP_2piRp(mask_ninf));% RLA
0084             %
0085             % power absorption rate
0086             %
0087             sP_2piRp_h = [sP_2piRp(1),(sP_2piRp(1:end-1) + sP_2piRp(2:end))/2,sP_2piRp(end)];% RLA
0088             ss_h = [ray.ss(1),(ray.ss(1:end-1) + ray.ss(2:end))/2,ray.ss(end)];
0089             sdP_2piRp_ds = -diff(sP_2piRp_h)./diff(ss_h);% RLA
0090             %
0091             % absorption coefficient
0092             %
0093             salpha = sdP_2piRp_ds./sP_2piRp;% RLA
0094             sphi_norm = sqrt(abs(ray.sphi_xyz(1,:)).^2 + abs(ray.sphi_xyz(2,:)).^2 + abs(ray.sphi_xyz(3,:)).^2);
0095             salphaphi = salpha.*sphi_norm;% RLA
0096             %
0097             waves{iw}.rays{ib}.sP_2piRp{itn} = sP_2piRp;% RLA
0098             waves{iw}.rays{ib}.stau{itn} = stau;% RLA
0099             waves{iw}.rays{ib}.sdP_2piRp_ds{itn} = sdP_2piRp_ds;% RLA
0100             waves{iw}.rays{ib}.salphaphi{itn} = salphaphi;% RLA
0101             %
0102             if isfield(dke_out,'zP_0_2piRp_mod_coll'),
0103                 %
0104                 zP_end_2piRp_coll = max([0,zP_0_2piRp_mod_coll{ibb}(end) - Prf_end_2piRp/2]);%to account for all QL absorbed power (RLA + NRCA)
0105                 %
0106                 zS_fit = [ray.ss(1),zS{ibb},ray.ss(end)];
0107                 zP_0_2piRp_coll_fit = [ray.P0_2piRp,zP_0_2piRp_mod_coll{ibb},zP_end_2piRp_coll];% RLA + NRCA
0108                 %
0109                 % power along the ray
0110                 %
0111                 sP_2piRp_coll = interp1(zS_fit,zP_0_2piRp_coll_fit,waves{iw}.rays{ib}.ss,'pchip');% RLA + NRCA
0112                 %
0113                 % optical depth
0114                 %
0115                 mask_ninf = sP_2piRp_coll > 0;
0116                 stau_coll = Inf(1,ns);
0117                 stau_coll(mask_ninf) = log(ray.P0_2piRp./sP_2piRp_coll(mask_ninf));%RLA + NRCA
0118                 %
0119                 % power absorption rate
0120                 %
0121                 sP_2piRp_coll_h = [sP_2piRp(1),(sP_2piRp_coll(1:end-1) + sP_2piRp_coll(2:end))/2,sP_2piRp_coll(end)];% RLA + NRCA
0122                 ss_h = [ray.ss(1),(ray.ss(1:end-1) + ray.ss(2:end))/2,ray.ss(end)];
0123                 sdP_2piRp_coll_ds = -diff(sP_2piRp_coll_h)./diff(ss_h);% RLA + NRCA
0124                 %
0125                 % absorption coefficient
0126                 %
0127                 salpha_coll = sdP_2piRp_coll_ds./sP_2piRp_coll;% RLA + NRCA
0128                 sphi_norm = sqrt(abs(ray.sphi_xyz(1,:)).^2 + abs(ray.sphi_xyz(2,:)).^2 + abs(ray.sphi_xyz(3,:)).^2);
0129                 salphaphi_coll = salpha_coll.*sphi_norm;% RLA + NRCA
0130                 %
0131                 waves{iw}.rays{ib}.sP_2piRp_coll_luke{itn} = sP_2piRp_coll;% RLA + NRCA
0132                 waves{iw}.rays{ib}.stau_coll_luke{itn} = stau_coll;% RLA + NRCA
0133                 waves{iw}.rays{ib}.sdP_2piRp_coll_luke_ds{itn} = sdP_2piRp_coll_ds;% RLA + NRCA
0134                 waves{iw}.rays{ib}.salphaphi_coll_luke{itn} = salphaphi_coll;% RLA + NRCA
0135                 %
0136             end 
0137         end
0138     end
0139 end
0140 %
0141 for iw = 1:nw,
0142     waves{iw}.simul = dke_out.simul;
0143 end
0144 %