alphaphimex_jd

C/C++ source

0001 /*================================================================================================
0002     alphaphimex_jd.c    MEX file implementation of : 
0003     alphaphi_jd.c,  which is a c subroutine for the calculation of 
0004     the fully relativistic absorption coefficient 
0005     for arbitrary RF waves interacting with electrons
0006 
0007 => xalphaphi = alphaphimex_jd(xa,xb,xbeta,xNpar,xNperp,xepolp,xepolm,xepolz,n,rel_opt,ns,pperpmax);
0008 
0009     INPUTS : 
0010 
0011         REQUIRED
0012         
0013         - xa    : omega_pe^2/omega_rf^2 [1,nx]
0014         - xb    : omega_ce/omega_rf [1,nx]
0015         - xbeta : sqrt(kTe/mc2) [1,nx]
0016         - xnpar : parallel index of refraction [1,nx]
0017         - xnperp: perpendicular index of refraction [1,nx]
0018         - xepolp: polarization vector component in rotating direction p [1,nx]
0019         - xepolm: polarization vector component in rotating direction m [1,nx]
0020         - xepolz: polarization vector component in parallel direction z [1,nx]
0021         - n     : harmonic number [1,1] 
0022         
0023         OPTIONAL
0024            
0025         - rel_opt : option for non-relativistic (0) or relativistic (1) calculations (default: 0)
0026         - ns      : number of pperp integration steps [1,1] (default: 1000)
0027         - pperpmax: maximum value for pperp integration [1,1] (default: 10)
0028   Note: pperpmax valid only where xNpar >= 1; elsewhere, pperpmax chosen by ellipse max position.
0029 
0030     OUTPUTS :
0031 
0032         - xalphaphi: alpha/(omega_rf/clum/phi) [1,nx]
0033 
0034 by Joan Decker (CEA-DRFC) <joan.decker@cea.fr> and Yves Peysson (CEA-DRFC)  <yves.peysson@cea.fr> 
0035 
0036 ================================================================================================*/
0037 
0038 /*#include <math.h>*/
0039 #include "mex.h"
0040 #include <complex.h>
0041 #include <math.h>
0042 
0043 /* subroutines list */
0044 
0045 double alphaphi_jd();
0046 
0047 /* Input Arguments */
0048 
0049 #define XA_IN           prhs[0]
0050 #define XB_IN           prhs[1]
0051 #define XBETA_IN        prhs[2]
0052 #define XNPAR_IN        prhs[3]
0053 #define XNPERP_IN       prhs[4]
0054 #define XEPOLP_R_IN     prhs[5]
0055 #define XEPOLP_I_IN     prhs[6]
0056 #define XEPOLM_R_IN     prhs[7]
0057 #define XEPOLM_I_IN     prhs[8]
0058 #define XEPOLZ_R_IN     prhs[9]
0059 #define XEPOLZ_I_IN     prhs[10]
0060 #define N_IN            prhs[11]
0061 #define REL_OPT_IN      prhs[12]
0062 #define NS_IN           prhs[13]
0063 #define PPERPMAX_IN     prhs[14]
0064 #define PN_IN           prhs[15]
0065 #define MHU_IN          prhs[16]
0066 #define F0_IN           prhs[17]
0067 #define PSI_IN          prhs[18]
0068 
0069 /* Output Arguments */
0070 
0071 #define XALPHAPHI_OUT   plhs[0]
0072 
0073 void mexFunction( int nlhs, mxArray *plhs[], 
0074             int nrhs, const mxArray *prhs[] )
0075      
0076 {     
0077     
0078     int n,rel_opt,ns;
0079     double pperpmax;
0080     double *xa,*xb,*xbeta,*xnpar,*xnperp;
0081     double *xepolp_r,*xepolp_i,*xepolm_r,*xepolm_i,*xepolz_r,*xepolz_i;
0082     double *xalphaphi;
0083     
0084     double complex xepolp,xepolm,xepolz;
0085     
0086     int flag_f0 =  1;
0087     
0088     int nx,m;
0089     
0090     double *pn,*mhu,*Xf0,*PSI;
0091     int npn,nmhu;
0092     
0093     /* Check out the arguments */
0094 
0095     if (nlhs != 1) {
0096         mexErrMsgTxt("Wrong number of output arguments for alphaphi_fr_jd");
0097     }
0098     
0099     if (nrhs < 12) {
0100         mexErrMsgTxt("Not enough input arguments calling alphaphi_fr_jd");
0101     }
0102  
0103     if (nrhs > 19)  {
0104         mexErrMsgTxt("Too many input arguments calling alphaphi_fr_jd");
0105     }
0106     
0107     if ((nrhs > 15) && (nrhs != 19))  {
0108         mexErrMsgTxt("Wrong number of input arguments calling alphaphi_fr_jd");
0109     }
0110       
0111     if (nrhs < 13) {
0112         rel_opt = 0;
0113     } else {
0114         rel_opt = *mxGetPr(REL_OPT_IN);
0115     }    
0116     
0117     if (nrhs < 14) {
0118         ns = 1000; 
0119     } else {
0120         ns = *mxGetPr(NS_IN);
0121     }    
0122     
0123     if (nrhs < 15) {
0124         pperpmax = 10;
0125     } else {
0126         pperpmax = *mxGetPr(PPERPMAX_IN);
0127     }
0128     
0129     if (nrhs != 19) flag_f0 = 0;/* No wave damping on a non-maxwellian distribution */
0130 
0131     
0132     /* Get Input Arguments */
0133 
0134     xa = mxGetPr(XA_IN);
0135     xb = mxGetPr(XB_IN);
0136     xbeta = mxGetPr(XBETA_IN);
0137     xnpar = mxGetPr(XNPAR_IN);
0138     xnperp = mxGetPr(XNPERP_IN);
0139     xepolp_r = mxGetPr(XEPOLP_R_IN);
0140     xepolp_i = mxGetPr(XEPOLP_I_IN);
0141     xepolm_r = mxGetPr(XEPOLM_R_IN);
0142     xepolm_i = mxGetPr(XEPOLM_I_IN);
0143     xepolz_r = mxGetPr(XEPOLZ_R_IN);
0144     xepolz_i = mxGetPr(XEPOLZ_I_IN);
0145     
0146     n = *mxGetPr(N_IN);
0147     
0148     /* Create Output Arguments */
0149     
0150     nx = mxGetN(XA_IN);
0151     XALPHAPHI_OUT = mxCreateDoubleMatrix(1,nx,mxREAL);
0152     xalphaphi = mxGetPr(XALPHAPHI_OUT);  
0153     
0154     if (flag_f0 == 1) {/* non-Maxwellian electron distribution */
0155 
0156         pn = mxGetPr(PN_IN);
0157         mhu = mxGetPr(MHU_IN);
0158         Xf0 = mxGetPr(F0_IN);
0159         PSI = mxGetPr(PSI_IN);
0160         
0161         npn = fmax(mxGetM(PN_IN),mxGetN(PN_IN));
0162         nmhu = fmax(mxGetM(MHU_IN),mxGetN(MHU_IN));
0163         
0164     } else {/* Maxwellian electron distribution */
0165           
0166         Xf0 = mxCalloc(1,sizeof(double));          
0167         PSI = mxCalloc(1,sizeof(double));
0168         
0169         pn = mxCalloc(1,sizeof(double));
0170         npn = 1;
0171         
0172         mhu = mxCalloc(1,sizeof(double));
0173         nmhu = 1;  
0174     }
0175     
0176     /* Local test 
0177     
0178     xalphaphi[0] = 0.1;
0179     xalphaphi[1] = 0.2;
0180     
0181     mexPrintf("xa[0],xa[1] : %g,%g\n",xa[0],xa[1]);
0182     mexPrintf("xb[0],xb[1] : %g,%g\n",xb[0],xb[1]);
0183     mexPrintf("xbeta[0],xbeta[1] : %g,%g\n",xbeta[0],xbeta[1]);
0184     mexPrintf("xnpar[0],xnpar[1] : %g,%g\n",xnpar[0],xnpar[1]);
0185     mexPrintf("xnperp[0],xnperp[1] : %g,%g\n",xnperp[0],xnperp[1]);
0186     
0187     mexPrintf("n,rel_opt,ns,pperpmax : %i,%i,%i,%g\n",n,rel_opt,ns,pperpmax);*/
0188 
0189     for (m = 0; m < nx; m++) {
0190         
0191         xepolp = xepolp_r[m] + I*xepolp_i[m];
0192         xepolm = xepolm_r[m] + I*xepolm_i[m];
0193         xepolz = xepolz_r[m] + I*xepolz_i[m];
0194         
0195         /*mexPrintf("xepolp : %g + I*%g\n",creal(xepolp),cimag(xepolp));
0196         mexPrintf("xepolm : %g + I*%g\n",creal(xepolm),cimag(xepolm));
0197         mexPrintf("xepolz : %g + I*%g\n",creal(xepolz),cimag(xepolz));*/
0198     
0199         xalphaphi[m] = alphaphi_jd(xepolp,xepolm,xepolz,xa[m],xb[m],xbeta[m],xnpar[m],xnperp[m],n,rel_opt,ns,pperpmax,npn,pn,nmhu,mhu,Xf0,PSI);
0200     }
0201     
0202     
0203     if (flag_f0 == 0) {/* Maxwellian electron distribution */
0204           
0205         mxFree(Xf0);
0206         mxFree(PSI);
0207         mxFree(pn);
0208         mxFree(mhu); 
0209     }  
0210 }
0211 
0212 
0213