alphaphi_fr_jd

   ALPHAPHI_FR_JD 

   This function calculates the fully relativistic electron absorption coefficient,
   normalized to omega_rf/clum/phi, where omega_rf is the wave angular frequency
   (in rad/s), clum is the speed of light (in m/s) and phi is the normalized
   energy flow density. The result is an adimentional quantity.

   INPUTS : 

           - xa: omega_pe^2/omega_rf^2 [1,nx]
           - xb: omega_ce/omega_rf [1,nx]
           - xbeta: sqrt(kTe/mc2) [1,nx]
           - xNpar: parallel index of refraction [1,nx]
           - xNperp: perpendicular index of refraction [1,nx]
           - xepol_pmz: polarization vector in rotating coordinates [3,nx]
           - n: harmonic number [1,1] 
           - rel_opt: option for non-relativistic (0) or relativistic (1) calculations (default: 0)
           - ns: number of pperp integration steps [1,1] (default: 1000)
           - pperpmax: maximum value for pperp integration [1,1] (default: 10)
               Note: pperpmax valid only where xNpar >= 1; 
                     elsewhere, pperpmax chosen by ellipse max position.
           - mex: (1) use or (0) not of the mex version if it exists

   OUTPUTS :

           - xalphaphi: alpha/(omega_rf/clum/phi) [1,nx]

   by J. Decker 11/17/05

0001 function xalphaphi = alphaphi_fr_jd(xa,xb,xbeta,xNpar,xNperp,xepol_pmz,n,rel_opt,ns,pperpmax,mex);
0002 %
0003 %   ALPHAPHI_FR_JD
0004 %
0005 %   This function calculates the fully relativistic electron absorption coefficient,
0006 %   normalized to omega_rf/clum/phi, where omega_rf is the wave angular frequency
0007 %   (in rad/s), clum is the speed of light (in m/s) and phi is the normalized
0008 %   energy flow density. The result is an adimentional quantity.
0009 %
0010 %   INPUTS :
0011 %
0012 %           - xa: omega_pe^2/omega_rf^2 [1,nx]
0013 %           - xb: omega_ce/omega_rf [1,nx]
0014 %           - xbeta: sqrt(kTe/mc2) [1,nx]
0015 %           - xNpar: parallel index of refraction [1,nx]
0016 %           - xNperp: perpendicular index of refraction [1,nx]
0017 %           - xepol_pmz: polarization vector in rotating coordinates [3,nx]
0018 %           - n: harmonic number [1,1]
0019 %           - rel_opt: option for non-relativistic (0) or relativistic (1) calculations (default: 0)
0020 %           - ns: number of pperp integration steps [1,1] (default: 1000)
0021 %           - pperpmax: maximum value for pperp integration [1,1] (default: 10)
0022 %               Note: pperpmax valid only where xNpar >= 1;
0023 %                     elsewhere, pperpmax chosen by ellipse max position.
0024 %           - mex: (1) use or (0) not of the mex version if it exists
0025 %
0026 %   OUTPUTS :
0027 %
0028 %           - xalphaphi: alpha/(omega_rf/clum/phi) [1,nx]
0029 %
0030 %   by J. Decker 11/17/05
0031 %
0032 if nargin < 11,
0033     mex = 1;
0034 end
0035 if nargin < 10,
0036     pperpmax = 10;
0037 end
0038 if nargin < 9,
0039     ns = 1000;
0040 end
0041 if nargin < 8,
0042     rel_opt = 0;
0043 end
0044 if nargin < 7,
0045     error('Not enough arguments')
0046 end
0047 %
0048 xa = xa(:).';
0049 xb = xb(:).';
0050 xbeta = xbeta(:).';
0051 xNpar = xNpar(:).';
0052 xNperp = xNperp(:).';
0053 %
0054 if size(xepol_pmz,2) == 3,
0055     if size(xepol_pmz,1) == 3,
0056         warning('Dimension of xepol_pmz is ambiguous');
0057     else
0058         xepol_pmz = xepol_pmz.';
0059     end
0060 end
0061 %
0062 nx = max([length(xa),length(xb),length(xbeta),length(xNpar),length(xNperp),size(xepol_pmz,2)]);
0063 %
0064 if length(xa) ~= nx,
0065     if length(xa) == 1,
0066         xa = xa*ones(1,nx);
0067     else
0068         error('Argument dimension mismatch')
0069     end
0070 end
0071 %
0072 if length(xb) ~= nx,
0073     if length(xb) == 1,
0074         xb = xb*ones(1,nx);
0075     else
0076         error('Argument dimension mismatch')
0077     end
0078 end
0079 %
0080 if length(xbeta) ~= nx,
0081     if length(xbeta) == 1,
0082         xbeta = xbeta*ones(1,nx);
0083     else
0084         error('Argument dimension mismatch')
0085     end
0086 end
0087 %
0088 if length(xNpar) ~= nx,
0089     if length(xNpar) == 1,
0090         xNpar = xNpar*ones(1,nx);
0091     else
0092         error('Argument dimension mismatch')
0093     end
0094 end
0095 %
0096 if length(xNperp) ~= nx,
0097     if length(xNperp) == 1,
0098         xNperp = xNperp*ones(1,nx);
0099     else
0100         error('Argument dimension mismatch')
0101     end
0102 end
0103 %
0104 if size(xepol_pmz,2) ~= nx,
0105     if size(xepol_pmz,2) == 1,
0106         xepol_pmz = xepol_pmz*ones(1,nx);
0107     else
0108         error('Argument dimension mismatch')
0109     end
0110 end
0111 %
0112 if size(xepol_pmz,1) ~= 3,
0113     error('Argument dimension mismatch')
0114 end
0115 %
0116 % USE OF MEX FILE IF IT EXISTS
0117 %
0118 if mex == 1 & exist('alphaphimex_jd'),
0119     xalphaphi = alphaphimex_jd(xa,xb,xbeta,xNpar,real(xNperp),...
0120         real(xepol_pmz(1,:)),imag(xepol_pmz(1,:)),...
0121         real(xepol_pmz(2,:)),imag(xepol_pmz(2,:)),...
0122         real(xepol_pmz(3,:)),imag(xepol_pmz(3,:)),...
0123         n,double(rel_opt),ns,pperpmax); 
0124 %        xepol_pmz(1,:),xepol_pmz(2,:),xepol_pmz(3,:),...
0125     return
0126 end
0127 %
0128 % MATLAB IMPLEMENTATION
0129 %
0130 xmask = ~isnan(xa) & ~isnan(xb) & ~isnan(xbeta) & ~isnan(xNpar) & ~isnan(xNperp) & ~isnan(sum(xepol_pmz,1));
0131 %
0132 xyn = xb*n;
0133 xslperp = xNperp.*xbeta./xb;
0134 %
0135 Ntest = ceil(min(2*pi*ns/pperpmax./xslperp));
0136 if Ntest < 50,
0137     warning(['Some positions have as little as ',num2str(Ntest),' integration steps per Bessel period'])
0138 end
0139 %
0140 xmask_ellipse = xmask & (rel_opt == 1) & (xNpar.^2 < 1) & (n > 0) & (xyn > sqrt(1 - xNpar.^2));
0141 xmask_parabola = xmask & (rel_opt == 1) & (xNpar.^2 == 1) & (n > 0);
0142 xmask_hyperbola = xmask & (rel_opt == 1) & (xNpar.^2 > 1);
0143 xmask_straight = xmask & (rel_opt == 0) & ((1 - xyn).^2 < xNpar.^2);
0144 %
0145 xpperpmax = pperpmax*ones(1,nx); % maximum pperp value in the integration
0146 xpperpmax_ellipse = sqrt(xyn.^2./(1 - xNpar.^2) - 1)./xbeta; % maximum value of the ellipse
0147 xpperpmax(xmask_ellipse) = xpperpmax_ellipse(xmask_ellipse);
0148 xpperpmax(xpperpmax > pperpmax) = pperpmax; % pperp limited to pperp <= pperpmax
0149 %
0150 xInt = zeros(1,nx);
0151 xInt(~xmask) = NaN;
0152 %
0153 if nx*ns <= 1e7,%no loop needed on the position
0154     %
0155     xInt(xmask_ellipse) = integral_fr_jd(xpperpmax(1,xmask_ellipse),xNpar(1,xmask_ellipse),xyn(1,xmask_ellipse),xbeta(1,xmask_ellipse),xepol_pmz(1,xmask_ellipse),xepol_pmz(2,xmask_ellipse),xepol_pmz(3,xmask_ellipse),xslperp(1,xmask_ellipse),n,ns,-1)...
0156                         + integral_fr_jd(xpperpmax(1,xmask_ellipse),xNpar(1,xmask_ellipse),xyn(1,xmask_ellipse),xbeta(1,xmask_ellipse),xepol_pmz(1,xmask_ellipse),xepol_pmz(2,xmask_ellipse),xepol_pmz(3,xmask_ellipse),xslperp(1,xmask_ellipse),n,ns,+1);
0157     %
0158     xInt(xmask_parabola) = integral_fr_jd(xpperpmax(1,xmask_parabola),xNpar(1,xmask_parabola),xyn(1,xmask_parabola),xbeta(1,xmask_parabola),xepol_pmz(1,xmask_parabola),xepol_pmz(2,xmask_parabola),xepol_pmz(3,xmask_parabola),xslperp(1,xmask_parabola),n,ns,0);
0159     %
0160     xInt(xmask_hyperbola) = integral_fr_jd(xpperpmax(1,xmask_hyperbola),xNpar(1,xmask_hyperbola),xyn(1,xmask_hyperbola),xbeta(1,xmask_hyperbola),xepol_pmz(1,xmask_hyperbola),xepol_pmz(2,xmask_hyperbola),xepol_pmz(3,xmask_hyperbola),xslperp(1,xmask_hyperbola),n,ns,-1);
0161     %
0162     xInt(xmask_straight) = integral_fr_jd(xpperpmax(1,xmask_straight),xNpar(1,xmask_straight),xyn(1,xmask_straight),xbeta(1,xmask_straight),xepol_pmz(1,xmask_straight),xepol_pmz(2,xmask_straight),xepol_pmz(3,xmask_straight),xslperp(1,xmask_straight),n,ns,NaN);
0163     %
0164 else % loop needed on the position
0165     %
0166     %h = waitbar(0,'Calculation of the absorption coefficient');
0167     %
0168     for ix = 1:nx,
0169             %
0170         if xmask_ellipse(ix) == 1,
0171             %
0172             xInt(ix) = integral_fr_jd(xpperpmax(ix),xNpar(ix),xyn(ix),xbeta(ix),xepol_pmz(1,ix),xepol_pmz(2,ix),xepol_pmz(3,ix),xslperp(ix),n,ns,-1)...
0173                      + integral_fr_jd(xpperpmax(ix),xNpar(ix),xyn(ix),xbeta(ix),xepol_pmz(1,ix),xepol_pmz(2,ix),xepol_pmz(3,ix),xslperp(ix),n,ns,+1);
0174             %
0175         elseif xmask_parabola(ix) == 1,
0176             %
0177             xInt(ix) = integral_fr_jd(xpperpmax(ix),xNpar(ix),xyn(ix),xbeta(ix),xepol_pmz(1,ix),xepol_pmz(2,ix),xepol_pmz(3,ix),xslperp(ix),n,ns,0);
0178             %
0179         elseif xmask_hyperbola(ix) == 1,
0180             %
0181             xInt(ix) = integral_fr_jd(xpperpmax(ix),xNpar(ix),xyn(ix),xbeta(ix),xepol_pmz(1,ix),xepol_pmz(2,ix),xepol_pmz(3,ix),xslperp(ix),n,ns,-1);
0182             %
0183         elseif xmask_straight(ix) == 1,
0184             %
0185             xInt(ix) = integral_fr_jd(xpperpmax(ix),xNpar(ix),xyn(ix),xbeta(ix),xepol_pmz(1,ix),xepol_pmz(2,ix),xepol_pmz(3,ix),xslperp(ix),n,ns,NaN);
0186             %
0187         end  
0188         %
0189         %waitbar(ix/nx,h);
0190         %
0191     end
0192     %
0193     %close(h);
0194     %
0195 end
0196 %
0197 if rel_opt == 1,
0198     xR = sqrt(pi/2)*xbeta./besselk(2,1./xbeta.^2,1);
0199 else
0200     xR = 1;
0201 end
0202 %
0203 xalphaphi = sqrt(pi/2)*xa.*xR.*xInt./xbeta;
0204 %
0205 %   Subfunction for the perpendicular momentum space integral
0206 %
0207 function xInt = integral_fr_jd(xpperpmax,xNpar,xyn,xbeta,xepol_p,xepol_m,xepol_z,xslperp,n,ns,res_opt);
0208     %
0209     sxpperp_S = linspace(0,1,ns+1).'*xpperpmax;   
0210     sxdpperp = diff(sxpperp_S);
0211     sxpperp = (sxpperp_S(1:ns,:) + sxpperp_S(2:ns+1,:))/2;
0212     %
0213     clear sxpperp_S
0214     %
0215     xsigma = sign(xNpar);
0216     %
0217     sxNpar = repmat(xNpar,[ns,1]);
0218     sxsigma = repmat(xsigma,[ns,1]);
0219     sxyn = repmat(xyn,[ns,1]);
0220     sxbeta = repmat(xbeta,[ns,1]);
0221     %
0222     if res_opt < 0,%branch m
0223         sxppar_res = (sxNpar.*sxyn - sxsigma.*sqrt(sxyn.^2 - (1 - sxNpar.^2).*(1 + sxbeta.^2.*sxpperp.^2)))./(sxbeta.*(1 - sxNpar.^2));
0224     elseif res_opt > 0,%branch p
0225         sxppar_res = (sxNpar.*sxyn + sxsigma.*sqrt(sxyn.^2 - (1 - sxNpar.^2).*(1 + sxbeta.^2.*sxpperp.^2)))./(sxbeta.*(1 - sxNpar.^2));
0226     elseif res_opt == 0,%branch 0
0227         sxppar_res = sxsigma.*(1 - sxyn.^2 + sxpperp.^2.*sxbeta.^2)./(2*sxbeta.*sxyn);
0228     else % non-relativistic calculation
0229         sxppar_res = (1 - sxyn)./sxbeta./sxNpar;
0230     end
0231     %
0232     clear sxsigma sxyn 
0233     %
0234     sxepol_p =  repmat(xepol_p,[ns,1]);
0235     sxepol_m =  repmat(xepol_m,[ns,1]);
0236     sxepol_z =  repmat(xepol_z,[ns,1]);
0237     %
0238     sxzeta = sxpperp.*repmat(xslperp,[ns,1]);
0239     %
0240     sxtheta = (sxepol_p.*besselj(n+1,sxzeta) + sxepol_m.*besselj(n-1,sxzeta))/sqrt(2) ...
0241         + sxepol_z.*besselj(n,sxzeta).*sxppar_res./sxpperp;
0242     %
0243     if ~isnan(res_opt),% relativistic case
0244         sxgamma = sqrt(1 + (sxppar_res.^2 + sxpperp.^2).*sxbeta.^2);
0245         %
0246         sxfint = sxpperp.^3.*abs(sxtheta).^2.*exp(-(sxgamma - 1)./sxbeta.^2)./abs(sxgamma.*sxNpar - sxppar_res.*sxbeta);
0247         %
0248     else
0249         %
0250         sxfint = sxpperp.^3.*abs(sxtheta).^2.*exp(-(sxppar_res.^2 + sxpperp.^2)/2)./abs(sxNpar);
0251         %
0252     end        
0253     %
0254     xInt = sum(sxfint.*sxdpperp);
0255 %