alphaphi_jd

C/C++ source

0001 #include "mex.h"
0002 #include <complex.h>
0003 #include <math.h>
0004 
0005 #include <bessel_IJK_jd.h>
0006 
0007 /* Format definition */
0008 
0009 
0010 double integral_fr_jd();
0011 
0012 #define PI 3.141592653589793
0013 
0014 double alphaphi_jd(epolp,epolm,epolz,a,b,beta,npar,nperp,n,rel_opt,ns,pperpmax,npn,pn,nmhu,mhu,Xf0,PSI)
0015 double complex epolp,epolm,epolz;
0016 double a,b,beta,npar,nperp,pperpmax;
0017 int n,rel_opt,ns;
0018 int npn,nmhu;
0019 double *pn,*mhu,*Xf0,*PSI;
0020 {   /* Local test */
0021      
0022     /*mexPrintf("a : %g\n",a);
0023     mexPrintf("b : %g\n",b);
0024     mexPrintf("beta : %g\n",beta);
0025     mexPrintf("npar : %g\n",npar);
0026     mexPrintf("nperp : %g\n",nperp);
0027     mexPrintf("epolp : %g + i*%g\n",creal(epolp),cimag(epolp));
0028     mexPrintf("epolm : %g + i*%g\n",creal(epolm),cimag(epolm));
0029     mexPrintf("epolz : %g + i*%g\n",creal(epolz),cimag(epolz));
0030     mexPrintf("n,rel_opt,ns,pperpmax : %i,%i,%i,%g\n",n,rel_opt,ns,pperpmax);*/
0031 
0032     /*return(a);*/
0033     
0034     double absint,pperpmax_ellipse,r;
0035     double yn,slperp,npar2,beta2,ntest;
0036     
0037     yn = b*n;
0038     slperp = nperp*beta/b;
0039     npar2 = npar*npar;
0040     beta2 = beta*beta;
0041     
0042     /* Tests for pperp grid accuracy */
0043     
0044     /* mexPrintf("ns,pperpmax,slperp : %i,%g,%g\n",ns,pperpmax,slperp);*/
0045     ntest = 2*PI*ns/pperpmax/slperp;
0046     
0047     /* mexPrintf("ntest : %g\n",ntest);*/
0048     
0049     if (ntest < 50) {
0050         mexPrintf("WARNING : Some positions have as little as %g integration steps per Bessel period\n",floor(ntest));
0051     }
0052     
0053     /* Non-Maxwellian parameter */
0054     
0055     if ((npn > 1) && (nmhu > 1)) {/*Wave aborption is calculated with a non-Maxwellian distribution function */
0056         
0057        /* mexPrintf("----------------------------------------\n");
0058         mexPrintf("Wave aborption may be calculated with a non-Maxwellian distribution function.\n"); 
0059         mexPrintf("npn, pn(1),pn(end) : %d,%g,%g\n",npn,pn[0],pn[npn-1]);
0060         mexPrintf("nmhu, mhu(1),mhu(end) : %d,%g,%g\n",nmhu,mhu[0],mhu[nmhu-1]);
0061         */
0062   
0063     }
0064       
0065     /* Resonance curve type */
0066     
0067     absint = 0;
0068     
0069     if (rel_opt == 1) {/* relativistic case */
0070         if (npar2 > 1) {/* hyperbola */
0071             absint = integral_fr_jd(epolp,epolm,epolz,pperpmax,npar,yn,beta,slperp,n,ns,-1); 
0072         } else {        
0073         if (npar2 == 1 && n > 0) {/* parabola */
0074             absint = integral_fr_jd(epolp,epolm,epolz,pperpmax,npar,yn,beta,slperp,n,ns,0); 
0075         } else {        
0076         if (npar2 < 1 && n > 0 && yn > sqrt(1 - npar2)) {/* ellipse */
0077             
0078             /* pperpmax redefined for ellipse */
0079             
0080             pperpmax_ellipse = sqrt(yn*yn/(1 - npar2) - 1)/beta;
0081             if (pperpmax_ellipse > pperpmax) {
0082                 pperpmax_ellipse = pperpmax;
0083             }
0084             
0085             absint = integral_fr_jd(epolp,epolm,epolz,pperpmax_ellipse,npar,yn,beta,slperp,n,ns,1); 
0086             absint = absint + integral_fr_jd(epolp,epolm,epolz,pperpmax_ellipse,npar,yn,beta,slperp,n,ns,-1); 
0087         }}}
0088     } else {/* non-relativistic case - straight line */
0089         if ((1 - yn)*(1 - yn) < npar2) {
0090             absint = integral_fr_jd(epolp,epolm,epolz,pperpmax,npar,yn,beta,slperp,n,ns,2); 
0091         }
0092     }
0093 
0094     
0095     if (rel_opt == 1) {
0096         if (beta2 < 0.002) { /* asymptotic expression */
0097             r = 1 - 15*beta2/8 + 345*beta2*beta2/128;
0098         } else {
0099             r = sqrt(PI/2)*beta/besskn(2,1/beta2)/exp(1/beta2);
0100         }
0101         
0102         /*mexPrintf("besskn(2,1/beta2)) : %g\n",besskn(2,1/beta2));
0103         mexPrintf("exp(1/beta2) : %g\n",exp(1/beta2));
0104         mexPrintf("r : %g\n",r);*/
0105         
0106     } else {
0107         r = 1;
0108     }
0109     
0110     /*mexPrintf("beta2,absint : %g,%g\n",beta2,absint);*/
0111     
0112     return(sqrt(PI/2)*a*r*absint/beta);
0113 }
0114 
0115 double integral_fr_jd(epolp,epolm,epolz,pperpmax,npar,yn,beta,slperp,n,ns,resopt) 
0116 double complex epolp,epolm,epolz;
0117 double pperpmax,npar,yn,beta,slperp;
0118 int n,ns,resopt;
0119 {
0120     /* simple method extracted from MATLAB version */
0121     
0122     double absint,fint,dpperp,pperp,pperp2,pparres,pparres2,zeta,theta2,gamma;
0123     double complex theta;
0124     double yn2,npar2,beta2;
0125     int m,sigma;
0126     
0127     /* mexPrintf("pperpmax,npar,yn,beta,epolp,epolm,epolz,slperp,n,ns,resopt : %g,%g,%g,%g,%g,%g,%g,%g,%i,%i,%i\n",pperpmax,npar,yn,beta,epolp,epolm,epolz,slperp,n,ns,resopt); */
0128     
0129     dpperp = pperpmax/ns;
0130     
0131     /* mexPrintf("ns,pperpmax,dpperp : %i,%g,%g\n",ns,pperpmax,dpperp);*/
0132     
0133     pperp = dpperp/2;
0134     yn2 = yn*yn;
0135     npar2 = npar*npar;
0136     beta2 = beta*beta;
0137     
0138     absint = 0;
0139     
0140     for (m = 0; m < ns; m++) {
0141         
0142         pperp2 = pperp*pperp;
0143         
0144         /* resonance condition gives ppar*/
0145         
0146         sigma = 0;
0147         if (npar > 0) {
0148             sigma = 1;
0149         } else {
0150         if (npar < 0) { 
0151             sigma = -1;
0152         }}        
0153         
0154         if (resopt == -1) {/* branch m */
0155             pparres = (npar*yn - sigma*sqrt(yn2 - (1 - npar2)*(1 + beta2*pperp2)))/(beta*(1 - npar2));            
0156         } else {
0157         if (resopt == 1) {/* branch p */
0158             pparres = (npar*yn + sigma*sqrt(yn2 - (1 - npar2)*(1 + beta2*pperp2)))/(beta*(1 - npar2));
0159         } else {
0160         if (resopt == 0) {/* branch 0 */
0161             pparres = sigma*(1 - yn2 + pperp2*beta2)/(2*beta*yn);
0162         } else { /* non-relativistic calculation */
0163             pparres = (1 - yn)/beta/npar;
0164         }}}
0165 
0166         /* mexPrintf("m,pparres,npar,npar2,sigma,yn,yn2,beta,beta2,pperp2 : %i,%g,%g,%g,%i,%g,%g,%g,%g,%g\n",m,pparres,npar,npar2,sigma,yn,yn2,beta,beta2,pperp2); */
0167 
0168         pparres2 = pparres*pparres;
0169         
0170         /* polarization term */
0171         
0172         zeta = pperp*slperp;
0173         theta = (epolp*bessjn(n+1,zeta) + epolm*bessjn(n-1,zeta))/sqrt(2) + epolz*bessjn(n,zeta)*pparres/pperp;
0174         theta2 = cabs(theta)*cabs(theta);
0175         
0176         /* integrand */
0177         
0178         if (resopt*resopt <= 1) { /* relativistic case */
0179             gamma = sqrt(1 + (pparres2 + pperp2)*beta2);
0180             fint = pperp2*pperp*theta2*exp(-(gamma - 1)/beta2)/fabs(gamma*npar - pparres*beta);
0181             
0182             /* mexPrintf("m,pparres,pperp,theta2,gamma,fint : %i,%g,%g,%g,%g,%g\n",m,pparres,pperp,theta2,gamma,fint); */
0183         } else {
0184             fint = pperp2*pperp*theta2*exp(-(pparres2 + pperp2)/2)/fabs(npar);
0185             
0186              /* mexPrintf("m,pparres,pperp,theta2,npar,fint : %i,%g,%g,%g,%g,%g\n",m,pparres,pperp,theta2,npar,fint); */
0187         }
0188 
0189         absint = absint + fint*dpperp;
0190         
0191         /* mexPrintf("m,pparres,pperp,fint,absint : %i,%g,%g,%g,%g\n",m,pparres,pperp,fint,absint); */
0192 
0193         pperp = pperp + dpperp;
0194                 
0195     }
0196         
0197     return(absint);    
0198 }    
0199