disp_hot_jd

 Hot plasma dispersion solver R2D2 

 by J. Decker and A.K. Ram

0001 function [smask_ht,sNperp_ht,sepol_pmz_ht,sPflow_pmz_ht,sPabs_at,...
0002         smask_ft,sNperp_ft,sepol_pmz_ft] = ...
0003     disp_hot_jd(ss,sB,sne,sTe,sNpar,omega,sNperp,rel_opt,tol_ht,nd_ht,guess_ES,nn_nr);
0004 %
0005 % Hot plasma dispersion solver R2D2
0006 %
0007 % by J. Decker and A.K. Ram
0008 %
0009 nnpar = 50;%number of points to continuously raise Npar when Npar = 1.
0010 %
0011 if nargin < 12
0012     nn_nr = 50;
0013 end
0014 if nargin < 11
0015     guess_ES = 10;
0016 end
0017 if nargin < 10
0018     nd_ht = 11;
0019 end
0020 if nargin < 9
0021     tol_ht = 0.2;
0022 end
0023 if nargin < 8
0024     rel_opt = 0;%non-relativistic calculations
0025 end
0026 %
0027 if nargin < 7
0028     sNperp = NaN*sNpar;%electrostatic root
0029 end
0030 %
0031 ns = length(ss);
0032 %
0033 [qe,me,mp,mn,e0,mu0,re,mc2,clum] = pc_dke_yp;
0034 %
0035 if sum(isnan(sNperp)) == ns% look for electrostatic root
0036     %
0037     % Cut-Offs location in cold plasma
0038     %
0039     swpe = sqrt(sne*qe^2/e0/me);
0040     swce = qe*sB/me;
0041     %
0042     iOco_app = min(find(omega < swpe));
0043     %
0044     if sNpar(iOco_app) < 1,
0045         swR = sqrt(swce.^2/4 + swpe.^2./(1-sNpar.^2)) + swce/2;
0046         swL = sqrt(swce.^2/4 + swpe.^2./(1-sNpar.^2)) - swce/2; 
0047         %
0048         iRco_app = min(find(omega < swR));
0049         iLco_app = min(find(omega < swL));
0050         iCmin_app = max([iRco_app,iLco_app,iOco_app]);
0051         %mask_npar = [];
0052     else % in this case, we continously search roots going from Npar = Npar_0 to the final value
0053         %sNpar_guess = 1/sqrt(1 + omega/swce(iOco_app));%Npar_0
0054         %mask_npar = [iOco_app]*ones(1,nnpar);
0055         iCmin_app = iOco_app;
0056     end    
0057     %
0058     % ES approximation - guess for hot plasma calculations
0059     %
0060     gn = 1 - round(omega/swce(iCmin_app))*swce(iCmin_app)/omega;
0061     if gn > 0,
0062         guess_Nperp = 50;%LFR approach
0063     else
0064         guess_Nperp = -50*gn;%HFR approach
0065     end
0066     %
0067     fac_root = 1.1;
0068     fac_exp = -1.5;
0069     fac_lim = 20;
0070     GNperp = NaN;
0071     %
0072     while (isnan(GNperp) | real(GNperp) < 1) & fac_root < fac_lim
0073         [GNperp] = disp_ebwES_jd(sB(iCmin_app),sTe(iCmin_app),sne(iCmin_app),omega,NaN,guess_Nperp*fac_root);
0074         fac_root = fac_root^fac_exp;
0075     end
0076     %
0077     % Determines ESW propagation range
0078     %
0079     mask_in = iCmin_app:ns;
0080     %
0081 else% look for electromagnetic root - closest from input Nperp
0082     %
0083     GNperp = sNperp(1);
0084     %
0085     mask_in = 1:ns;
0086 end
0087 %
0088 ns_in = length(mask_in);
0089 ss_in = ss(mask_in);
0090 sNpar_in = sNpar(mask_in);
0091 sB_in = sB(mask_in);
0092 sTe_in = sTe(mask_in);
0093 sne_in = sne(mask_in);
0094 %
0095 knr = 1-rel_opt;
0096 nrt_nr = 1;
0097 nrt_fr = 1;
0098 kflag = 1;
0099 %
0100 [ss_out_ft,sMNperp_ft,sM3eps_pmz_ft,ss_out_ht,sMNperp_ht,sM3eps_pmz_ht,sM3Pflow_xyz_ht,sMPabs_at,sM3Pflow_pmz_ht] = ...
0101     disp_jd(ss_in,sB_in,sTe_in,sne_in,sNpar_in,omega,knr,nrt_nr,nrt_fr,kflag,GNperp,nn_nr,tol_ht,nd_ht);
0102 ns_out_ft = length(ss_out_ft);
0103 ns_out_ht = length(ss_out_ht);
0104 if ns_out_ht ~= ns_in
0105     diff_ns = num2str(100 - ns_out_ht*100/ns_in);
0106     disp(['Warning: disp. calc. have been truncated by ',diff_ns(1:2),' %'])
0107 end
0108 %
0109 for is = 1:ns_out_ht,
0110     smask_ht(1,is) = find(abs(ss_out_ht(is) - ss) == min(abs(ss_out_ht(is) - ss))); 
0111 end
0112 for is = 1:ns_out_ft,
0113     smask_ft(1,is) = find(abs(ss_out_ft(is) - ss) == min(abs(ss_out_ft(is) - ss))); 
0114 end
0115 sNperp_ht = sMNperp_ht';%since nrt=1
0116 sepol_pmz_ht = reshape(sM3eps_pmz_ht,[ns_out_ht,3])';
0117 %sPflow_xyz_ht = reshape(sM3Pflow_xyz_ht,[ns_out_ht,3])';
0118 sPflow_pmz_ht = reshape(sM3Pflow_pmz_ht,[ns_out_ht,3])';
0119 sPabs_at = sMPabs_at';
0120 sNperp_ft = sMNperp_ft';
0121 sepol_pmz_ft = reshape(sM3eps_pmz_ft,[ns_out_ft,3])';
0122 %
0123