raytracing_dieltensor_yp

C/C++ source

0001 #include <complex.h>
0002 #include <math.h>
0003 #include "mex.h"
0004 
0005 
0006 
0007 
0008 
0009 struct raymagnetic {
0010     double *rn;
0011     double *drndY;
0012     double *d2rndY2;
0013     double *xn;
0014     double *dxndY;
0015     double *d2xndY2;
0016     double *yn;
0017     double *dyndY;
0018     double *d2yndY2;
0019     double *Brhon;
0020     double *dBrhondY;
0021     double *d2BrhondY2;
0022     double *Bsn;
0023     double *dBsndY;
0024     double *d2BsndY2;
0025     double *Bzn;
0026     double *dBzndY;
0027     double *d2BzndY2;
0028     double *B;
0029     double *dBdY;
0030     double *d2BdY2;
0031     double *salpha;
0032     double *dsalphadY;
0033     double *d2salphadY2;
0034     double *calpha;
0035     double *dcalphadY;
0036     double *d2calphadY2;
0037     double *gradrhon;
0038     double *dgradrhondY;
0039     double *d2gradrhondY2;
0040     double *rhorip;
0041     double *drhoripdY;
0042     double *d2rhoripdY2;
0043     double *Upsilonn;
0044     double *dUpsilonndY;
0045     double *d2UpsilonndY2;
0046 };
0047 
0048 
0049 
0050 struct susceptibilitytensor {
0051     double complex *Xxx;
0052     double complex *Xxy;
0053     double complex *Xxz;
0054     double complex *Xyy;
0055     double complex *Xyz;
0056     double complex *Xzz;
0057     double complex *dXxxdNperp;
0058     double complex *dXxydNperp;
0059     double complex *dXxzdNperp;
0060     double complex *dXyydNperp;
0061     double complex *dXyzdNperp;
0062     double complex *dXzzdNperp;
0063     double complex *dXxxdNpar;
0064     double complex *dXxydNpar;
0065     double complex *dXxzdNpar;
0066     double complex *dXyydNpar;
0067     double complex *dXyzdNpar;
0068     double complex *dXzzdNpar;
0069     double complex *dXxxdY;
0070     double complex *dXxydY;
0071     double complex *dXxzdY;
0072     double complex *dXyydY;
0073     double complex *dXyzdY;
0074     double complex *dXzzdY; 
0075 };
0076 
0077 typedef struct raymagnetic raymagnetic_format;
0078 typedef struct susceptibilitytensor susceptibilitytensor_format;
0079 
0080 
0081 #define MC 8
0082 #define MAXCYCLHARM 101
0083 #define PI 3.14159265358979
0084 
0085 void buildXstruc();
0086 void vacuumtensor();
0087 void coldtensor();
0088 void calcY();
0089 void calcY1();
0090 void Zdisp();
0091 void ZFriedConte();
0092 void hottensor();
0093 void Gamma();
0094 void modBessel();
0095 void imagBessel();
0096 void dawson();
0097 void smith();
0098 void weakrelativistictensor();
0099 void Shka();
0100 void ImShka();
0101 void Fdnest();
0102 void ImFdnest();
0103 int sgn();
0104 
0105 
0106 void vacuumtensor(X)
0107 susceptibilitytensor_format X;
0108 { 
0109     int i;
0110     
0111     X.Xxx[0] = 0.0;
0112     X.Xxy[0] = 0.0;
0113     X.Xxz[0] = 0.0;
0114     X.Xyy[0] = 0.0;
0115     X.Xyz[0] = 0.0;
0116     X.Xzz[0] = 0.0;
0117 
0118     for (i = 0;i < MC;i++) {
0119          X.dXxxdY[i] = 0.0;
0120          X.dXxydY[i] = 0.0;
0121          X.dXxzdY[i] = 0.0;
0122          X.dXyydY[i] = 0.0;
0123          X.dXyzdY[i] = 0.0;
0124          X.dXzzdY[i] = 0.0;
0125     }
0126 } 
0127 
0128 void coldtensor(X,Nperp,dNperpdY,ns,ncyclharm,raymagnetic,Npar,dNpardY,wcn,dwcndY,wpn2,dwpn2dY,betath,dbetathdY)
0129 susceptibilitytensor_format X;
0130 double complex Nperp,dNperpdY[];
0131 int ns;
0132 int ncyclharm;
0133 raymagnetic_format raymagnetic;
0134 double wcn[],wpn2[],betath[];
0135 double dwcndY[][MC],dwpn2dY[][MC],dbetathdY[][MC];
0136 double Npar,dNpardY[];
0137 
0138 { 
0139     double complex Xxx[ns+1],Xxy[ns+1],Xxz[ns+1],Xyy[ns+1],Xyz[ns+1],Xzz[ns+1];
0140     double complex dXxxdNperp[ns+1],dXxydNperp[ns+1],dXxzdNperp[ns+1],dXyydNperp[ns+1],dXyzdNperp[ns+1],dXzzdNperp[ns+1];
0141     double complex dXxxdNpar[ns+1],dXxydNpar[ns+1],dXxzdNpar[ns+1],dXyydNpar[ns+1],dXyzdNpar[ns+1],dXzzdNpar[ns+1];
0142     double complex dXxxdbetath[ns+1],dXxydbetath[ns+1],dXxzdbetath[ns+1],dXyydbetath[ns+1],dXyzdbetath[ns+1],dXzzdbetath[ns+1];
0143     double complex dXxxdwcn[ns+1],dXxydwcn[ns+1],dXxzdwcn[ns+1],dXyydwcn[ns+1],dXyzdwcn[ns+1],dXzzdwcn[ns+1];
0144     double complex dXxxdwpn2[ns+1],dXxydwpn2[ns+1],dXxzdwpn2[ns+1],dXyydwpn2[ns+1],dXyzdwpn2[ns+1],dXzzdwpn2[ns+1];
0145     
0146     int is,i;
0147 
0148     calcY1(Nperp,dNperpdY,ns,ncyclharm,Npar,dNpardY,wcn,dwcndY,wpn2,dwpn2dY,betath,dbetathdY);/* For problems on some LINUX systems with gcc */
0149     
0150     for (is = 0;is < ns + 1;is++) {
0151         
0152 /* Dielectric tensor */
0153         
0154         Xxx[is] = - wpn2[is]/(1.0 - wcn[is]*wcn[is]);
0155         Xxy[is] = I*wpn2[is]*wcn[is]/(1.0 - wcn[is]*wcn[is]);
0156         Xxz[is] = 0.0;
0157         Xyy[is] = - wpn2[is]/(1.0 - wcn[is]*wcn[is]);
0158         Xyz[is] = 0.0;
0159         Xzz[is] = - wpn2[is];
0160 
0161 /* derivatives of the dielectric tensor */
0162         
0163         dXxxdNperp[is] = 0.0;
0164         dXxydNperp[is] = 0.0;
0165         dXxzdNperp[is] = 0.0;
0166         dXyydNperp[is] = 0.0;
0167         dXyzdNperp[is] = 0.0;
0168         dXzzdNperp[is] = 0.0;
0169 
0170         dXxxdNpar[is] = 0.0;
0171         dXxydNpar[is] = 0.0;
0172         dXxzdNpar[is] = 0.0;
0173         dXyydNpar[is] = 0.0;
0174         dXyzdNpar[is] = 0.0;
0175         dXzzdNpar[is] = 0.0;
0176 
0177         dXxxdbetath[is] = 0.0;
0178         dXxydbetath[is] = 0.0;
0179         dXxzdbetath[is] = 0.0;
0180         dXyydbetath[is] = 0.0;
0181         dXyzdbetath[is] = 0.0;
0182         dXzzdbetath[is] = 0.0;
0183 
0184         dXxxdwcn[is] = -2.0*wpn2[is]*wcn[is]/(1.0 - wcn[is]*wcn[is])/(1.0 - wcn[is]*wcn[is]);
0185         dXxydwcn[is] = I*(2.0*wpn2[is]*wcn[is]*wcn[is]/(1 - wcn[is]*wcn[is])/(1.0 - wcn[is]*wcn[is]) + wpn2[is]/(1.0 - wcn[is]*wcn[is]));
0186         dXxzdwcn[is] = 0.0;
0187         dXyydwcn[is] = -2.0*wpn2[is]*wcn[is]/(1.0 - wcn[is]*wcn[is])/(1.0 - wcn[is]*wcn[is]);
0188         dXyzdwcn[is] = 0.0;
0189         dXzzdwcn[is] = 0.0;    
0190 
0191         dXxxdwpn2[is] = -1.0/(1.0 - wcn[is]*wcn[is]);
0192         dXxydwpn2[is] = I*wcn[is]/(1.0 - wcn[is]*wcn[is]);
0193         dXxzdwpn2[is] = 0.0;
0194         dXyydwpn2[is] = -1.0/(1.0 - wcn[is]*wcn[is]);
0195         dXyzdwpn2[is] = 0.0;
0196         dXzzdwpn2[is] = -1.0;    
0197         
0198     }
0199      
0200     buildXstruc(X,dNperpdY,Xxx,Xxy,Xxz,Xyy,Xyz,Xzz,dXxxdNperp,dXxydNperp,dXxzdNperp,dXyydNperp,dXyzdNperp,dXzzdNperp,dXxxdNpar,dXxydNpar,dXxzdNpar,dXyydNpar,dXyzdNpar,dXzzdNpar,dXxxdbetath,dXxydbetath,dXxzdbetath,dXyydbetath,dXyzdbetath,dXzzdbetath,dXxxdwcn,dXxydwcn,dXxzdwcn,dXyydwcn,dXyzdwcn,dXzzdwcn,dXxxdwpn2,dXxydwpn2,dXxzdwpn2,dXyydwpn2,dXyzdwpn2,dXzzdwpn2,ns,dNpardY,dwcndY,dwpn2dY,dbetathdY);
0201 } 
0202 
0203 
0204 
0205 void calcY1(Nperp,dNperpdY,ns,ncyclharm,Npar,dNpardY,wcn,dwcndY,wpn2,dwpn2dY,betath,dbetathdY)
0206 double complex Nperp,dNperpdY[MC];
0207 int ns;
0208 int ncyclharm;
0209 double Npar,dNpardY[MC];
0210 double wcn[],dwcndY[][MC],wpn2[],dwpn2dY[][MC],betath[],dbetathdY[][MC];
0211 { 
0212 
0213 }
0214 
0215 
0216 void hottensor(X,Nperp,dNperpdY,ns,ncyclharm,raymagnetic,Npar,dNpardY,wcn,dwcndY,wpn2,dwpn2dY,betath,dbetathdY)
0217 susceptibilitytensor_format X;
0218 double complex Nperp,dNperpdY[];
0219 int ns;
0220 int ncyclharm;
0221 raymagnetic_format raymagnetic;
0222 double wcn[],wpn2[],betath[];
0223 double dwcndY[][MC],dwpn2dY[][MC],dbetathdY[][MC];
0224 double Npar,dNpardY[];
0225 
0226 { 
0227     double complex Xxx[ns+1],Xxy[ns+1],Xxz[ns+1],Xyy[ns+1],Xyz[ns+1],Xzz[ns+1];
0228     double complex dXxxdNperp[ns+1],dXxydNperp[ns+1],dXxzdNperp[ns+1],dXyydNperp[ns+1],dXyzdNperp[ns+1],dXzzdNperp[ns+1];
0229     double complex dXxxdNpar[ns+1],dXxydNpar[ns+1],dXxzdNpar[ns+1],dXyydNpar[ns+1],dXyzdNpar[ns+1],dXzzdNpar[ns+1];
0230     double complex dXxxdbetath[ns+1],dXxydbetath[ns+1],dXxzdbetath[ns+1],dXyydbetath[ns+1],dXyzdbetath[ns+1],dXzzdbetath[ns+1];
0231     double complex dXxxdwcn[ns+1],dXxydwcn[ns+1],dXxzdwcn[ns+1],dXyydwcn[ns+1],dXyzdwcn[ns+1],dXzzdwcn[ns+1];
0232     double complex dXxxdwpn2[ns+1],dXxydwpn2[ns+1],dXxzdwpn2[ns+1],dXyydwpn2[ns+1],dXyzdwpn2[ns+1],dXzzdwpn2[ns+1];
0233     
0234     double Yxx[ns+1][2*ncyclharm+1],Yxy[ns+1][2*ncyclharm+1],Yxz[ns+1][2*ncyclharm+1],Yyy[ns+1][2*ncyclharm+1],Yyz[ns+1][2*ncyclharm+1],Yzz[ns+1][2*ncyclharm+1];
0235     double dYxxdNperp[ns+1][2*ncyclharm+1],dYxydNperp[ns+1][2*ncyclharm+1],dYxzdNperp[ns+1][2*ncyclharm+1],dYyydNperp[ns+1][2*ncyclharm+1],dYyzdNperp[ns+1][2*ncyclharm+1],dYzzdNperp[ns+1][2*ncyclharm+1];
0236     double dYxxdNpar[ns+1][2*ncyclharm+1],dYxydNpar[ns+1][2*ncyclharm+1],dYxzdNpar[ns+1][2*ncyclharm+1],dYyydNpar[ns+1][2*ncyclharm+1],dYyzdNpar[ns+1][2*ncyclharm+1],dYzzdNpar[ns+1][2*ncyclharm+1];
0237     double dYxxdbetath[ns+1][2*ncyclharm+1],dYxydbetath[ns+1][2*ncyclharm+1],dYxzdbetath[ns+1][2*ncyclharm+1],dYyydbetath[ns+1][2*ncyclharm+1],dYyzdbetath[ns+1][2*ncyclharm+1],dYzzdbetath[ns+1][2*ncyclharm+1];
0238     double dYxxdwcn[ns+1][2*ncyclharm+1],dYxydwcn[ns+1][2*ncyclharm+1],dYxzdwcn[ns+1][2*ncyclharm+1],dYyydwcn[ns+1][2*ncyclharm+1],dYyzdwcn[ns+1][2*ncyclharm+1],dYzzdwcn[ns+1][2*ncyclharm+1];
0239     double dYxxdwpn2[ns+1][2*ncyclharm+1],dYxydwpn2[ns+1][2*ncyclharm+1],dYxzdwpn2[ns+1][2*ncyclharm+1],dYyydwpn2[ns+1][2*ncyclharm+1],dYyzdwpn2[ns+1][2*ncyclharm+1],dYzzdwpn2[ns+1][2*ncyclharm+1];
0240     
0241     double zeta,dzetadNperp,dzetadNpar,dzetadbetath,dzetadwcn,dzetadwpn2; 
0242     
0243     int is,i;
0244 
0245     calcY(Nperp,dNperpdY,Yxx,Yxy,Yxz,Yyy,Yyz,Yzz,dYxxdNperp,dYxydNperp,dYxzdNperp,dYyydNperp,dYyzdNperp,dYzzdNperp,dYxxdNpar,dYxydNpar,dYxzdNpar,dYyydNpar,dYyzdNpar,dYzzdNpar,dYxxdbetath,dYxydbetath,dYxzdbetath,dYyydbetath,dYyzdbetath,dYzzdbetath,dYxxdwcn,dYxydwcn,dYxzdwcn,dYyydwcn,dYyzdwcn,dYzzdwcn,dYxxdwpn2,dYxydwpn2,dYxzdwpn2,dYyydwpn2,dYyzdwpn2,dYzzdwpn2,ns,ncyclharm,Npar,dNpardY,wcn,dwcndY,wpn2,dwpn2dY,betath,dbetathdY);
0246 
0247     for (is = 0;is < ns + 1;is++) {/* initialisation */           
0248         Xxx[is] = 0.0;
0249         Xxy[is] = 0.0;
0250         Xxz[is] = 0.0;
0251         Xyy[is] = 0.0;
0252         Xyz[is] = 0.0;
0253         Xzz[is] = 0.0;
0254 
0255         dXxxdNperp[is] = 0.0;
0256         dXxydNperp[is] = 0.0;
0257         dXxzdNperp[is] = 0.0;
0258         dXyydNperp[is] = 0.0;
0259         dXyzdNperp[is] = 0.0;
0260         dXzzdNperp[is] = 0.0;
0261 
0262         dXxxdNpar[is] = 0.0;
0263         dXxydNpar[is] = 0.0;
0264         dXxzdNpar[is] = 0.0;
0265         dXyydNpar[is] = 0.0;
0266         dXyzdNpar[is] = 0.0;
0267         dXzzdNpar[is] = 0.0;
0268 
0269         dXxxdbetath[is] = 0.0;
0270         dXxydbetath[is] = 0.0;
0271         dXxzdbetath[is] = 0.0;
0272         dXyydbetath[is] = 0.0;
0273         dXyzdbetath[is] = 0.0;
0274         dXzzdbetath[is] = 0.0;
0275 
0276         dXxxdwcn[is] = 0.0;
0277         dXxydwcn[is] = 0.0;
0278         dXxzdwcn[is] = 0.0;
0279         dXyydwcn[is] = 0.0;
0280         dXyzdwcn[is] = 0.0;
0281         dXzzdwcn[is] = 0.0;    
0282 
0283         dXxxdwpn2[is] = 0.0;
0284         dXxydwpn2[is] = 0.0;
0285         dXxzdwpn2[is] = 0.0;
0286         dXyydwpn2[is] = 0.0;
0287         dXyzdwpn2[is] = 0.0;
0288         dXzzdwpn2[is] = 0.0;    
0289     }
0290     
0291     for (is = 0;is < ns + 1;is++) {
0292         for (i = -ncyclharm;i <= ncyclharm; i++) {
0293             
0294             zeta = (1.0 - i*wcn[is])/fabs(Npar)*betath[is];
0295             dzetadNperp = 0.0;
0296             dzetadNpar = -sgn(Npar)*zeta/Npar;
0297             dzetadbetath = -zeta/betath[is];
0298             dzetadwcn = - i/fabs(Npar)*betath[is];
0299             dzetadwpn2 = 0.0;
0300             
0301 /* Dielectric tensor */
0302             
0303             Xxx[is] = Xxx[is] + wpn2[is]*zeta*Yxx[is][i];
0304             Xxy[is] = Xxy[is] + wpn2[is]*zeta*Yxy[is][i];
0305             Xxz[is] = Xxz[is] + wpn2[is]*zeta*Yxz[is][i];
0306             Xyy[is] = Xyy[is] + wpn2[is]*zeta*Yyy[is][i];
0307             Xyz[is] = Xyz[is] + wpn2[is]*zeta*Yyz[is][i];
0308             Xzz[is] = Xzz[is] + wpn2[is]*zeta*Yzz[is][i];
0309 
0310  /* derivatives of the dielectric tensor */
0311 
0312             dXxxdNperp[is] = dXxxdNperp[is] + wpn2[is]*dzetadNperp*Yxx[is][i] + wpn2[is]*zeta*dYxxdNperp[is][i];
0313             dXxydNperp[is] = dXxydNperp[is] + wpn2[is]*dzetadNperp*Yxy[is][i] + wpn2[is]*zeta*dYxydNperp[is][i];
0314             dXxzdNperp[is] = dXxzdNperp[is] + wpn2[is]*dzetadNperp*Yxz[is][i] + wpn2[is]*zeta*dYxzdNperp[is][i];
0315             dXyydNperp[is] = dXyydNperp[is] + wpn2[is]*dzetadNperp*Yyy[is][i] + wpn2[is]*zeta*dYyydNperp[is][i];
0316             dXyzdNperp[is] = dXyzdNperp[is] + wpn2[is]*dzetadNperp*Yyz[is][i] + wpn2[is]*zeta*dYyzdNperp[is][i];
0317             dXzzdNperp[is] = dXzzdNperp[is] + wpn2[is]*dzetadNperp*Yzz[is][i] + wpn2[is]*zeta*dYzzdNperp[is][i];
0318 
0319             dXxxdNpar[is] = dXxxdNpar[is] + wpn2[is]*dzetadNpar*Yxx[is][i] + wpn2[is]*zeta*dYxxdNpar[is][i];
0320             dXxydNpar[is] = dXxydNpar[is] + wpn2[is]*dzetadNpar*Yxy[is][i] + wpn2[is]*zeta*dYxydNpar[is][i];
0321             dXxzdNpar[is] = dXxzdNpar[is] + wpn2[is]*dzetadNpar*Yxz[is][i] + wpn2[is]*zeta*dYxzdNpar[is][i];
0322             dXyydNpar[is] = dXyydNpar[is] + wpn2[is]*dzetadNpar*Yyy[is][i] + wpn2[is]*zeta*dYyydNpar[is][i];
0323             dXyzdNpar[is] = dXyzdNpar[is] + wpn2[is]*dzetadNpar*Yyz[is][i] + wpn2[is]*zeta*dYyzdNpar[is][i];
0324             dXzzdNpar[is] = dXzzdNpar[is] + wpn2[is]*dzetadNpar*Yzz[is][i] + wpn2[is]*zeta*dYzzdNpar[is][i];
0325 
0326             dXxxdbetath[is] = dXxxdbetath[is] + wpn2[is]*dzetadbetath*Yxx[is][i] + wpn2[is]*zeta*dYxxdbetath[is][i];
0327             dXxydbetath[is] = dXxydbetath[is] + wpn2[is]*dzetadbetath*Yxy[is][i] + wpn2[is]*zeta*dYxydbetath[is][i];
0328             dXxzdbetath[is] = dXxzdbetath[is] + wpn2[is]*dzetadbetath*Yxz[is][i] + wpn2[is]*zeta*dYxzdbetath[is][i];
0329             dXyydbetath[is] = dXyydbetath[is] + wpn2[is]*dzetadbetath*Yyy[is][i] + wpn2[is]*zeta*dYyydbetath[is][i];
0330             dXyzdbetath[is] = dXyzdbetath[is] + wpn2[is]*dzetadbetath*Yyz[is][i] + wpn2[is]*zeta*dYyzdbetath[is][i];
0331             dXzzdbetath[is] = dXzzdbetath[is] + wpn2[is]*dzetadbetath*Yzz[is][i] + wpn2[is]*zeta*dYzzdbetath[is][i];
0332 
0333             dXxxdwcn[is] = dXxxdwcn[is] + wpn2[is]*dzetadwcn*Yxx[is][i] + wpn2[is]*zeta*dYxxdwcn[is][i];
0334             dXxydwcn[is] = dXxydwcn[is] + wpn2[is]*dzetadwcn*Yxy[is][i] + wpn2[is]*zeta*dYxydwcn[is][i];
0335             dXxzdwcn[is] = dXxzdwcn[is] + wpn2[is]*dzetadwcn*Yxz[is][i] + wpn2[is]*zeta*dYxzdwcn[is][i];
0336             dXyydwcn[is] = dXyydwcn[is] + wpn2[is]*dzetadwcn*Yyy[is][i] + wpn2[is]*zeta*dYyydwcn[is][i];
0337             dXyzdwcn[is] = dXyzdwcn[is] + wpn2[is]*dzetadwcn*Yyz[is][i] + wpn2[is]*zeta*dYyzdwcn[is][i];
0338             dXzzdwcn[is] = dXzzdwcn[is] + wpn2[is]*dzetadwcn*Yzz[is][i] + wpn2[is]*zeta*dYzzdwcn[is][i];    
0339 
0340             dXxxdwpn2[is] = dXxxdwpn2[is] + zeta*Yxx[is][i] + wpn2[is]*dzetadwpn2*Yxx[is][i] + wpn2[is]*zeta*dYxxdwpn2[is][i];
0341             dXxydwpn2[is] = dXxydwpn2[is] + zeta*Yxy[is][i] + wpn2[is]*dzetadwpn2*Yxy[is][i] + wpn2[is]*zeta*dYxydwpn2[is][i];
0342             dXxzdwpn2[is] = dXxzdwpn2[is] + zeta*Yxz[is][i] + wpn2[is]*dzetadwpn2*Yxz[is][i] + wpn2[is]*zeta*dYxzdwpn2[is][i];
0343             dXyydwpn2[is] = dXyydwpn2[is] + zeta*Yyy[is][i] + wpn2[is]*dzetadwpn2*Yyy[is][i] + wpn2[is]*zeta*dYyydwpn2[is][i];
0344             dXyzdwpn2[is] = dXyzdwpn2[is] + zeta*Yyz[is][i] + wpn2[is]*dzetadwpn2*Yyz[is][i] + wpn2[is]*zeta*dYyzdwpn2[is][i];
0345             dXzzdwpn2[is] = dXzzdwpn2[is] + zeta*Yzz[is][i] + wpn2[is]*dzetadwpn2*Yzz[is][i] + wpn2[is]*zeta*dYzzdwpn2[is][i];    
0346         }
0347     }
0348     
0349     buildXstruc(X,dNperpdY,Xxx,Xxy,Xxz,Xyy,Xyz,Xzz,dXxxdNperp,dXxydNperp,dXxzdNperp,dXyydNperp,dXyzdNperp,dXzzdNperp,dXxxdNpar,dXxydNpar,dXxzdNpar,dXyydNpar,dXyzdNpar,dXzzdNpar,dXxxdbetath,dXxydbetath,dXxzdbetath,dXyydbetath,dXyzdbetath,dXzzdbetath,dXxxdwcn,dXxydwcn,dXxzdwcn,dXyydwcn,dXyzdwcn,dXzzdwcn,dXxxdwpn2,dXxydwpn2,dXxzdwpn2,dXyydwpn2,dXyzdwpn2,dXzzdwpn2,ns,dNpardY,dwcndY,dwpn2dY,dbetathdY);
0350 } 
0351 
0352 
0353 void calcY(Nperp,dNperpdY,Yxx,Yxy,Yxz,Yyy,Yyz,Yzz,dYxxdNperp,dYxydNperp,dYxzdNperp,dYyydNperp,dYyzdNperp,dYzzdNperp,dYxxdNpar,dYxydNpar,dYxzdNpar,dYyydNpar,dYyzdNpar,dYzzdNpar,dYxxdbetath,dYxydbetath,dYxzdbetath,dYyydbetath,dYyzdbetath,dYzzdbetath,dYxxdwcn,dYxydwcn,dYxzdwcn,dYyydwcn,dYyzdwcn,dYzzdwcn,dYxxdwpn2,dYxydwpn2,dYxzdwpn2,dYyydwpn2,dYyzdwpn2,dYzzdwpn2,ns,ncyclharm,Npar,dNpardY,wcn,dwcndY,wpn2,dwpn2dY,betath,dbetathdY)
0354 double complex Nperp,dNperpdY[MC];
0355 double complex Yxx[][MAXCYCLHARM],Yxy[][MAXCYCLHARM],Yxz[][MAXCYCLHARM],Yyy[][MAXCYCLHARM],Yyz[][MAXCYCLHARM],Yzz[][MAXCYCLHARM];
0356 double complex dYxxdNperp[][MAXCYCLHARM],dYxydNperp[][MAXCYCLHARM],dYxzdNperp[][MAXCYCLHARM],dYyydNperp[][MAXCYCLHARM],dYyzdNperp[][MAXCYCLHARM],dYzzdNperp[][MAXCYCLHARM];
0357 double complex dYxxdNpar[][MAXCYCLHARM],dYxydNpar[][MAXCYCLHARM],dYxzdNpar[][MAXCYCLHARM],dYyydNpar[][MAXCYCLHARM],dYyzdNpar[][MAXCYCLHARM],dYzzdNpar[][MAXCYCLHARM];
0358 double complex dYxxdbetath[][MAXCYCLHARM],dYxydbetath[][MAXCYCLHARM],dYxzdbetath[][MAXCYCLHARM],dYyydbetath[][MAXCYCLHARM],dYyzdbetath[][MAXCYCLHARM],dYzzdbetath[][MAXCYCLHARM];
0359 double complex dYxxdwcn[][MAXCYCLHARM],dYxydwcn[][MAXCYCLHARM],dYxzdwcn[][MAXCYCLHARM],dYyydwcn[][MAXCYCLHARM],dYyzdwcn[][MAXCYCLHARM],dYzzdwcn[][MAXCYCLHARM];
0360 double complex dYxxdwpn2[][MAXCYCLHARM],dYxydwpn2[][MAXCYCLHARM],dYxzdwpn2[][MAXCYCLHARM],dYyydwpn2[][MAXCYCLHARM],dYyzdwpn2[][MAXCYCLHARM],dYzzdwpn2[][MAXCYCLHARM];
0361 int ns;
0362 int ncyclharm;
0363 double Npar,dNpardY[MC];
0364 double wcn[],dwcndY[][MC],wpn2[],dwpn2dY[][MC],betath[],dbetathdY[][MC];
0365 { 
0366     double zeta,dzetadNperp,dzetadNpar,dzetadbetath,dzetadwcn,dzetadwpn2;
0367     double complex lambda,dlambdadNperp,dlambdadNpar,dlambdadbetath,dlambdadwcn,dlambdadwpn2;
0368     double complex Z,Zp,Zpp,G,Gp,Gpp;
0369     double sigmas;
0370     
0371     int is,i;  
0372  
0373     if (2*ncyclharm+1 > MAXCYCLHARM) {mexPrintf("******** WARNING: the code is limited to 50 cylotron harmonics only ! ********");}
0374     
0375     for (is = 0;is < ns + 1;is++) {
0376         
0377         sigmas = sgn(-1.0+2.0*is);/* sign of the charge, electron is always first, the ions */
0378         
0379         for (i = -ncyclharm;i <= ncyclharm; i++) {
0380             
0381             zeta = (1.0 - i*wcn[is])/fabs(Npar)*betath[is];
0382             dzetadNperp = 0.0;
0383             dzetadNpar = -sgn(Npar)*zeta/Npar;
0384             dzetadbetath = -zeta/betath[is];
0385             dzetadwcn = - i/fabs(Npar)*betath[is];
0386             dzetadwpn2 = 0.0;
0387             
0388             lambda = Nperp*betath[is]/wcn[is];
0389             dlambdadNperp = betath[is]/wcn[is];
0390             dlambdadNpar = 0.0;
0391             dlambdadbetath = Nperp/wcn[is];
0392             dlambdadwcn = - lambda/wcn[is];
0393             dlambdadwpn2 = 0.0;
0394             
0395             Zdisp((double complex)zeta,&Z,&Zp,&Zpp);
0396             Gamma(lambda,&G,&Gp,&Gpp,(double)i);
0397             
0398 /* Dielectric tensor */
0399                        
0400             Yxx[is][i] = i*i*G*Z/lambda;
0401             Yxy[is][i] = I*sigmas*i*Gp*Z;
0402             Yxz[is][i] = -sgn(Npar)*i*G*Zp/sqrt(2.0*lambda);
0403             Yyy[is][i] = (i*i*G/lambda - 2.0*lambda*Gp)*Z;
0404             Yyz[is][i] = I*sigmas*sgn(Npar)*sqrt(lambda/2.0)*Gp*Zp;
0405             Yzz[is][i] = -G*zeta*Zp;
0406 
0407  /* derivatives of the dielectric tensor */
0408 
0409             dYxxdNperp[is][i] = i*i*Gp*dlambdadNperp*Z/lambda + i*i*G*Zp*dzetadNperp - i*i*G*Z*dlambdadNperp/lambda/lambda;
0410             dYxydNperp[is][i] = I*sigmas*i*Gpp*dlambdadNperp*Z + I*sigmas*i*Gp*Zp*dzetadNperp;
0411             dYxzdNperp[is][i] = -sgn(Npar)*i*Gp*dlambdadNperp*Zp/sqrt(2.0*lambda) - sgn(Npar)*i*G*Zpp*dzetadNperp/sqrt(2.0*lambda) + sgn(Npar)*i*G*Zp*dlambdadNperp/sqrt(2.0*lambda)/(2.0*lambda);
0412             dYyydNperp[is][i] = (i*i*Gp*dlambdadNperp/lambda - i*i*G*dlambdadNperp/lambda/lambda - 2.0*dlambdadNperp*Gp - 2.0*lambda*dlambdadNperp*Gpp)*Z + (i*i*G/lambda - 2.0*lambda*Gp)*Zp*dzetadNperp;
0413             dYyzdNperp[is][i] = I*sigmas*sgn(Npar)*dlambdadNperp*Gp*Zp/sqrt(lambda/2.0)/4.0 + I*sigmas*sgn(Npar)*sqrt(lambda/2.0)*Gpp*dlambdadNperp*Zp + I*sigmas*sgn(Npar)*sqrt(lambda/2.0)*Gp*Zpp*dzetadNperp;
0414             dYzzdNperp[is][i] = -Gp*dlambdadNperp*zeta*Zp - G*dzetadNperp*Zp - G*zeta*Zpp*dzetadNperp;
0415 
0416             dYxxdNpar[is][i] = i*i*Gp*dlambdadNpar*Z/lambda + i*i*G*Zp*dzetadNpar - i*i*G*Z*dlambdadNpar/lambda/lambda;
0417             dYxydNpar[is][i] = I*sigmas*i*Gpp*dlambdadNpar*Z + I*sigmas*i*Gp*Zp*dzetadNpar;
0418             dYxzdNpar[is][i] = -sgn(Npar)*i*Gp*dlambdadNpar*Zp/sqrt(2.0*lambda) - sgn(Npar)*i*G*Zpp*dzetadNpar/sqrt(2.0*lambda) + sgn(Npar)*i*G*Zp*dlambdadNpar/sqrt(2.0*lambda)/(2.0*lambda);
0419             dYyydNpar[is][i] = (i*i*Gp*dlambdadNpar/lambda - i*i*G*dlambdadNpar/lambda/lambda - 2.0*dlambdadNpar*Gp - 2.0*lambda*dlambdadNpar*Gpp)*Z + (i*i*G/lambda - 2.0*lambda*Gp)*Zp*dzetadNpar;
0420             dYyzdNpar[is][i] = I*sigmas*sgn(Npar)*dlambdadNpar*Gp*Zp/sqrt(lambda/2.0)/4.0 + I*sigmas*sgn(Npar)*sqrt(lambda/2.0)*Gpp*dlambdadNpar*Zp + I*sigmas*sgn(Npar)*sqrt(lambda/2.0)*Gp*Zpp*dzetadNpar;
0421             dYzzdNpar[is][i] = -Gp*dlambdadNpar*zeta*Zp - G*dzetadNpar*Zp - G*zeta*Zpp*dzetadNpar;
0422 
0423             dYxxdbetath[is][i] = i*i*Gp*dlambdadbetath*Z/lambda + i*i*G*Zp*dzetadbetath - i*i*G*Z*dlambdadbetath/lambda/lambda;
0424             dYxydbetath[is][i] = I*sigmas*i*Gpp*dlambdadbetath*Z + I*sigmas*i*Gp*Zp*dzetadbetath;
0425             dYxzdbetath[is][i] = -sgn(Npar)*i*Gp*dlambdadbetath*Zp/sqrt(2.0*lambda) - sgn(Npar)*i*G*Zpp*dzetadbetath/sqrt(2.0*lambda) + sgn(Npar)*i*G*Zp*dlambdadbetath/sqrt(2.0*lambda)/(2.0*lambda);
0426             dYyydbetath[is][i] = (i*i*Gp*dlambdadbetath/lambda - i*i*G*dlambdadbetath/lambda/lambda - 2.0*dlambdadbetath*Gp - 2.0*lambda*dlambdadbetath*Gpp)*Z + (i*i*G/lambda - 2.0*lambda*Gp)*Zp*dzetadbetath;
0427             dYyzdbetath[is][i] = I*sigmas*sgn(Npar)*dlambdadbetath*Gp*Zp/sqrt(lambda/2.0)/4.0 + I*sigmas*sgn(Npar)*sqrt(lambda/2.0)*Gpp*dlambdadbetath*Zp + I*sigmas*sgn(Npar)*sqrt(lambda/2.0)*Gp*Zpp*dzetadbetath;
0428             dYzzdbetath[is][i] = -Gp*dlambdadbetath*zeta*Zp - G*dzetadbetath*Zp - G*zeta*Zpp*dzetadbetath;
0429 
0430             dYxxdwcn[is][i] = i*i*Gp*dlambdadwcn*Z/lambda + i*i*G*Zp*dzetadwcn - i*i*G*Z*dlambdadwcn/lambda/lambda;
0431             dYxydwcn[is][i] = I*sigmas*i*Gpp*dlambdadwcn*Z + I*sigmas*i*Gp*Zp*dzetadwcn;
0432             dYxzdwcn[is][i] = -sgn(Npar)*i*Gp*dlambdadwcn*Zp/sqrt(2.0*lambda) - sgn(Npar)*i*G*Zpp*dzetadwcn/sqrt(2.0*lambda) + sgn(Npar)*i*G*Zp*dlambdadwcn/sqrt(2.0*lambda)/(2.0*lambda);
0433             dYyydwcn[is][i] = (i*i*Gp*dlambdadwcn/lambda - i*i*G*dlambdadwcn/lambda/lambda - 2.0*dlambdadwcn*Gp - 2.0*lambda*dlambdadwcn*Gpp)*Z + (i*i*G/lambda - 2.0*lambda*Gp)*Zp*dzetadwcn;
0434             dYyzdwcn[is][i] = I*sigmas*sgn(Npar)*dlambdadwcn*Gp*Zp/sqrt(lambda/2.0)/4.0 + I*sigmas*sgn(Npar)*sqrt(lambda/2.0)*Gpp*dlambdadwcn*Zp + I*sigmas*sgn(Npar)*sqrt(lambda/2.0)*Gp*Zpp*dzetadwcn;
0435             dYzzdwcn[is][i] = -Gp*dlambdadwcn*zeta*Zp - G*dzetadwcn*Zp - G*zeta*Zpp*dzetadwcn;    
0436 
0437             dYxxdwpn2[is][i] = i*i*Gp*dlambdadwpn2*Z/lambda + i*i*G*Zp*dzetadwpn2 - i*i*G*Z*dlambdadwpn2/lambda/lambda;
0438             dYxydwpn2[is][i] = I*sigmas*i*Gpp*dlambdadwpn2*Z + I*sigmas*i*Gp*Zp*dzetadwpn2;
0439             dYxzdwpn2[is][i] = -sgn(Npar)*i*Gp*dlambdadwpn2*Zp/sqrt(2.0*lambda) - sgn(Npar)*i*G*Zpp*dzetadwpn2/sqrt(2.0*lambda) + sgn(Npar)*i*G*Zp*dlambdadwpn2/sqrt(2.0*lambda)/(2.0*lambda);
0440             dYyydwpn2[is][i] = (i*i*Gp*dlambdadwpn2/lambda - i*i*G*dlambdadwpn2/lambda/lambda - 2.0*dlambdadwpn2*Gp - 2.0*lambda*dlambdadwpn2*Gpp)*Z + (i*i*G/lambda - 2.0*lambda*Gp)*Zp*dzetadwpn2;
0441             dYyzdwpn2[is][i] = I*sigmas*sgn(Npar)*dlambdadwpn2*Gp*Zp/sqrt(lambda/2.0)/4.0 + I*sigmas*sgn(Npar)*sqrt(lambda/2.0)*Gpp*dlambdadwpn2*Zp + I*sigmas*sgn(Npar)*sqrt(lambda/2.0)*Gp*Zpp*dzetadwpn2;
0442             dYzzdwpn2[is][i] = -Gp*dlambdadwpn2*zeta*Zp - G*dzetadwpn2*Zp - G*zeta*Zpp*dzetadwpn2;    
0443         }
0444     }   
0445 }
0446 
0447 void Gamma(z,g,gp,gpp,nhu)
0448 double complex z,*g,*gp,*gpp;
0449 int nhu;
0450 { 
0451     double complex mBn,mBnp,mBnpp;
0452     
0453     modBessel(z,&mBn,&mBnp,&mBnpp,nhu);
0454     
0455     g[0] = cexp(-z)*mBn;
0456     gp[0] = -g[0] + cexp(-z)*mBnp;
0457     gpp[0] = g[0] - cexp(-z)*(2.0*mBnp - mBnpp);
0458     
0459 }
0460 
0461 
0462 void Zdisp(z,Z,Zp,Zpp)
0463 double complex z;
0464 double complex *Z,*Zp,*Zpp;
0465 {  
0466 
0467     ZFriedConte(z,Z);
0468 
0469     Zp[0] = -2.0*(1.0 + z*Z[0]);
0470     Zpp[0] = 4.0*z + Z[0]*(4*z*z - 2.0);
0471 }
0472 
0473 void ZFriedConte(z,Z)
0474 double complex z;
0475 double complex *Z;
0476 {
0477     double daw,x,y;
0478     if (cimag(z) == 0.0) {/* pure real argument */
0479         x = creal(z);
0480         dawson(x,&daw);
0481         Z[0]=-2.0*daw+I*sqrt(PI)*exp(x*x);
0482     } else if (creal(z) == 0.0) {/* pure imaginary argument */
0483         y = cimag(z);
0484         Z[0] = I*sqrt(PI)*exp(y*y)*erfc(-y);
0485     } else {
0486         smith(z,Z); 
0487     }
0488 }
0489 
0490 
0491 void modBessel(z,mBn,mBnp,mBnpp,nhu)
0492 double complex z,*mBn,*mBnp,*mBnpp;
0493 double nhu;
0494 {   
0495     double complex iBn,iBnp,iBnpp;
0496     
0497     imagBessel(I*z,&iBn,&iBnp,&iBnpp,nhu);
0498 
0499     mBn[0] = cpow(-I,nhu)*iBn;
0500     mBnp[0] = -cpow(-I,nhu+1.0)*iBnp;
0501     mBnpp[0] = -cpow(-I,nhu)*iBnpp;
0502     
0503 /*     mexPrintf("modBesseln:%g,%g\n",mBn[0]); */
0504     
0505 }
0506 
0507 
0508 
0509 void imagBessel(z,iBn,iBnp,iBnpp,nhu)
0510 double complex z,*iBn,*iBnp,*iBnpp;
0511 double nhu;
0512 
0513 {   
0514     mxArray *rhs[3], *lhs[1];
0515     
0516     double  *pr,*pi;
0517     double complex Besseln,Besselnm1,Besselnm2;
0518     double *Bessel_r,*Bessel_i;
0519     
0520     rhs[0] = mxCreateString("J");
0521     rhs[1] = mxCreateDoubleMatrix(1, 1, mxREAL);
0522     rhs[2] = mxCreateDoubleMatrix(1, 1, mxCOMPLEX);
0523     
0524     pr = mxGetPr(rhs[2]);
0525     pr[0] = creal(z);
0526     pi = mxGetPi(rhs[2]);
0527     pi[0] = cimag(z);
0528     
0529     pr = mxGetPr(rhs[1]);
0530     pr[0] = nhu;
0531     mexCallMATLAB(1,lhs,3,rhs,"besselmx");/* call the MatLab built-in Bessel function: [w,ierr] = besselmx(double('J'),0,z) */
0532     Bessel_r = mxGetData(lhs[0]);
0533     Bessel_i = mxGetImagData(lhs[0]);
0534     Besseln = Bessel_r[0] + I*Bessel_i[0];
0535     
0536 /*     mexPrintf("Besseln:%g,%g\n",Besseln); */
0537     
0538     pr = mxGetPr(rhs[1]);
0539     pr[0] = nhu - 1.0;
0540     mexCallMATLAB(1,lhs,3,rhs,"besselmx");/* call the MatLab built-in Bessel function: [w,ierr] = besselmx(double('J'),0,z) */
0541     Bessel_r = mxGetData(lhs[0]);
0542     Bessel_i = mxGetImagData(lhs[0]);
0543     Besselnm1 = Bessel_r[0] + I*Bessel_i[0];
0544     
0545 /*     mexPrintf("Besseln:%g,%g\n",Besselnm1); */
0546     
0547     
0548     pr = mxGetPr(rhs[1]);
0549     pr[0] = nhu - 2.0;
0550     mexCallMATLAB(1,lhs,3,rhs,"besselmx");/* call the MatLab built-in Bessel function: [w,ierr] = besselmx(double('J'),0,z) */
0551     Bessel_r = mxGetData(lhs[0]);
0552     Bessel_i = mxGetImagData(lhs[0]);
0553     Besselnm2 = Bessel_r[0] + I*Bessel_i[0];
0554 
0555 /*     mexPrintf("Besseln:%g,%g\n",Besselnm2); */
0556   
0557     iBn[0] = Besseln;
0558     iBnp[0] = Besselnm1 - nhu*Besseln/z;
0559     iBnpp[0] = Besselnm2 + (2.0*nhu - 1.0)*Besselnm1/z + nhu*(nhu + 1.0)*Besseln/z/z;
0560     
0561     /* cleanup allocated memory */
0562     
0563     mxDestroyArray(rhs[0]);
0564     mxDestroyArray(rhs[1]);
0565     mxDestroyArray(rhs[2]);
0566     mxDestroyArray(lhs[0]);
0567 }
0568 
0569 
0570 
0571 void smith(z,Z) 
0572 double complex z, *Z;
0573 {
0574 
0575 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0576 % Z dispersion function
0577 % with complex argument zeta.
0578 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0579 % For real argument use also Dawson integral
0580 % with Tchebychev polynomials
0581 % Z=-2.0*dawson(zeta)+i*sqrt(pi)*exp(-zeta.*zeta);
0582 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0583 % For pure imaginary argument also use
0584 % Z(i*x)=i*sqrt(pi)*exp(x.*x).*erfc(-x)
0585 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0586 % Inspired by Gary R. Smith
0587 % Phys. Fluids 27 (6) june 1984 p. 1510 
0588 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0589 */
0590 
0591     double complex z2, z_2, som, ZZ;
0592     double rz,iz;
0593     double criterion, sigma; 
0594     int nn; 
0595 
0596 /* % epsilon=1e-10
0597 % b=4.74
0598 % nmax1=58
0599 % nmax2=25
0600 % nmax3=27
0601 
0602 % epsilon=1e-12 */
0603 
0604     static double b=5.51; /* boundary between formulas */
0605     static double nmax1=75; /* nmax for series expansion */
0606     static double nmax2=34; /* nmax for continued fraction */
0607     static double nmax3=27; /* nmax for asymptotic expansion */
0608     static double complex czero = 0.0;
0609     static double complex cone = 1.0;
0610     static double complex isqrtpi = I*1.772453850905516;
0611 
0612     z2 = z*z;
0613     som = czero;
0614     rz = fabs(creal(z));
0615     iz = fabs(cimag(z));
0616 
0617 /* Series expansion around zeta = 0 formula A10 */
0618     
0619     if ((rz <= b) & (iz <= 2.0-1.5*rz/b)) {
0620         for (nn = nmax1;nn > 0;nn--) {
0621             som = 1.0/nn*z2*1.0/(2.0*nn + 1.0 + som);
0622         }
0623         ZZ = cexp(-z2)*(isqrtpi - 2.0*z*(som + 1.0));
0624     }
0625 
0626 /* Continued fraction expansion formula A12 & A13 */
0627 
0628     if (((rz <= b) & (iz > 2.0-1.5*rz/b)) | ((rz > b) & (iz > 0.5))) {
0629 
0630         z_2 = cone/z2;
0631         for (nn = nmax2;nn > 0;nn--) {
0632             som = -0.5*nn*z_2/(som + 1.0);
0633         }
0634         ZZ = -1.0/(z*(som + 1.0));        
0635         if (cimag(z) < 0.0) {
0636             ZZ = ZZ + 2.0*isqrtpi*cexp(-z2);
0637         }
0638     }  
0639 
0640 /* Asymptotic expansion formula A15 */
0641     if ((rz > b) & (iz <= 0.5)) {
0642         z_2 = cone/z2;
0643         for (nn = nmax3;nn > 0;nn--) {
0644             som = (nn-0.5)*z_2*(som + 1.0);
0645         }
0646         ZZ = -(som + 1.0)/z;    
0647         criterion = 0.7853981633974483331/rz;
0648         sigma = (cimag(z) < criterion) + (cimag(z) < -criterion);
0649         if (sigma > 0.0) {
0650             ZZ = ZZ + sigma*isqrtpi*cexp(-z2);
0651         }
0652     }
0653 
0654     Z[0] = ZZ;
0655 } 
0656 
0657 
0658 
0659 void dawson(zeta,dwork) 
0660 double zeta;
0661 double *dwork;
0662 {
0663 
0664 /*     A FUNCTION TO EVALUATE DAWSON'S INTEGRAL,    
0665        Dawson(Zeta) = Integral ( t=0->Zeta ) Exp(t**2-Zeta**2) dt,
0666        = -0.5 * Real part of Fried & Conte Z function   
0667        USING A CHEBYCHEV POLYNOMIAL SERIES.  
0668 */
0669 
0670     static double half = 0.5; 
0671     static double one = 1.0; 
0672     static double two = 2.0; 
0673     static double ten = 10.0;
0674     
0675 /* nmax=    5         => 1e-2 precision
0676                10         => 1e-4
0677                15         => 1e-6
0678                20         => single
0679                50         => double 
0680 */
0681     static int nmax = 50;
0682     
0683     static double tcoef [51]={2.3770636862073187, -0.2214259431356961, -0.0232524179380925, 0.1327809189298409, -0.1325824891002065,
0684            0.0825446940740374,-0.0337827509774703, 0.0070874724137126, 0.0009814185303952,-0.0010948499244205,
0685            0.0001341521888280, 0.0001204181702739,-0.0000351287219323,-0.0000146145750095, 0.0000062641524139,
0686            0.0000022709548584,-0.0000010280349044,-0.0000004487400975, 0.0000001532631036, 0.0000001005306690,
0687           -0.0000000166366048,-0.0000000226490692,-0.0000000006425396, 0.0000000046258504, 0.0000000011681305,
0688           -0.0000000007317958,-0.0000000004540732, 0.0000000000395128, 0.0000000001175806, 0.0000000000278655,
0689           -0.0000000000190620,-0.0000000000134858,-0.0000000000001751, 0.0000000000033395, 0.0000000000013606,
0690           -0.0000000000003122,-0.0000000000004896,-0.0000000000001295, 0.0000000000000741, 0.0000000000000694,
0691            0.0000000000000122,-0.0000000000000135,-0.0000000000000101,-0.0000000000000013, 0.0000000000000022,
0692            0.0000000000000016, 0.0000000000000002,-0.0000000000000004,-0.0000000000000003, 0.0000000000000000,
0693            0.0000000000000001};
0694     double zeta2,t,twot,tn,tnm1,tnm2;
0695     int ind;
0696 
0697     zeta2 = zeta*zeta;
0698     t = (zeta2 - ten)/(zeta2 + ten);
0699     twot = two*t;
0700     tnm1 = one;
0701     tn = t;
0702     dwork[0] = half*tcoef[0]*tnm1 + tcoef[1]*tn;
0703     for (ind=2;ind <= nmax;ind++) {
0704         tnm2 = tnm1;
0705         tnm1 = tn;
0706         tn = twot*tnm1 - tnm2;
0707         dwork[0] += tcoef[ind]*tn;
0708     }
0709     dwork[0] *= zeta/(one + two*zeta2);  
0710 }
0711 
0712 void weakrelativistictensor(X,Nperp,dNperpdx,ns,ncyclharm,raymagnetic,Npar,dNpardx,wcn,dwcndx,wpn2,dwpn2dx,betath,dbetathdx)
0713 susceptibilitytensor_format X;
0714 double complex Nperp,dNperpdx[];
0715 int ns,ncyclharm;
0716 raymagnetic_format raymagnetic;
0717 double wcn[],wpn2[],betath[];
0718 double dwcndx[][MC],dwpn2dx[][MC],dbetathdx[][MC];
0719 double Npar,dNpardx[];
0720 
0721 { 
0722     
0723     double complex Xxx[ns+1],Xxy[ns+1],Xxz[ns+1],Xyy[ns+1],Xyz[ns+1],Xzz[ns+1];
0724     double complex dXxxdNperp[ns+1],dXxydNperp[ns+1],dXxzdNperp[ns+1],dXyydNperp[ns+1],dXyzdNperp[ns+1],dXzzdNperp[ns+1];
0725     double complex dXxxdNpar[ns+1],dXxydNpar[ns+1],dXxzdNpar[ns+1],dXyydNpar[ns+1],dXyzdNpar[ns+1],dXzzdNpar[ns+1];
0726     double complex dXxxdbetath[ns+1],dXxydbetath[ns+1],dXxzdbetath[ns+1],dXyydbetath[ns+1],dXyzdbetath[ns+1],dXzzdbetath[ns+1];
0727     double complex dXxxdwcn[ns+1],dXxydwcn[ns+1],dXxzdwcn[ns+1],dXyydwcn[ns+1],dXyzdwcn[ns+1],dXzzdwcn[ns+1];
0728     double complex dXxxdwpn2[ns+1],dXxydwpn2[ns+1],dXxzdwpn2[ns+1],dXyydwpn2[ns+1],dXyzdwpn2[ns+1],dXzzdwpn2[ns+1];
0729     
0730     int is,i;
0731 
0732     double npar2, te, B, XX, Y_1, Y_2, mu, halfmu, mu2, delta1, delta2, delta3, mudelta1, mudelta2, mudelta3;
0733     double halfnpar2, psi2, twopsi2, pipo1;
0734     double F50, F70, F90, dF50dz, dF70dz, dF90dz, dF30da, dF50da, dF70da, dF90da, dF110da;
0735     double I50, I70, I90, dI50dz, dI70dz, dI90dz, dI30da, dI50da, dI70da, dI90da, dI110da;
0736     double F51, F71, F91, dF51dz, dF71dz, dF91dz,         dF51da, dF71da, dF91da, dF111da;
0737     double I51, I71, I91, dI51dz, dI71dz, dI91dz,         dI51da, dI71da, dI91da, dI111da;
0738     double      F72, F92,         dF72dz, dF92dz,                 dF72da, dF92da, dF112da;
0739     double      I72, I92,         dI72dz, dI92dz,                 dI72da, dI92da, dI112da;
0740     double           F93,                 dF93dz,                         dF93da;
0741     double           I93,                 dI93dz,                         dI93da;
0742 
0743     double muF50, dmuF50dTe, dmuF50dnpar2, muI50, dmuI50dTe, dmuI50dnpar2;
0744     double dF70dTe, dF70dnpar2, dI70dTe, dI70dnpar2, dF90dTe, dF90dnpar2, dI90dTe, dI90dnpar2;
0745     double muF701, dmuF701dTe, dmuF701dnpar2, muI701, dmuI701dTe, dmuI701dnpar2;
0746     double F901, dF901dTe, dF901dnpar2, I901, dI901dTe, dI901dnpar2;
0747     double twopsi2muF702, dtwopsi2muF702dTe, dtwopsi2muF702dnpar2, twopsi2muI702, dtwopsi2muI702dTe, dtwopsi2muI702dnpar2;
0748     double twopsi2F902, dtwopsi2F902dTe, dtwopsi2F902dnpar2, twopsi2I902, dtwopsi2I902dTe, dtwopsi2I902dnpar2;
0749     double twopsi2F1102,dtwopsi2F1102dTe,dtwopsi2F1102dnpar2,twopsi2I1102,dtwopsi2I1102dTe,dtwopsi2I1102dnpar2;
0750     double muF51, dmuF51dY, dmuF51dTe, dmuF51dnpar2, muI51, dmuI51dY, dmuI51dTe, dmuI51dnpar2;
0751     double dF71dY, dF71dTe, dF71dnpar2, dI71dY, dI71dTe, dI71dnpar2;
0752     double dF91dY, dF91dTe, dF91dnpar2, dI91dY, dI91dTe, dI91dnpar2;
0753     double muF711, dmuF711dY, dmuF711dTe, dmuF711dnpar2, muI711, dmuI711dY, dmuI711dTe, dmuI711dnpar2;
0754     double F911, dF911dY, dF911dTe, dF911dnpar2, I911, dI911dY, dI911dTe, dI911dnpar2;
0755     double twopsi2F912,  dtwopsi2F912dY,  dtwopsi2F912dTe,  dtwopsi2F912dnpar2;
0756     double twopsi2I912,  dtwopsi2I912dY,  dtwopsi2I912dTe,  dtwopsi2I912dnpar2;
0757     double twopsi2F1112, dtwopsi2F1112dY, dtwopsi2F1112dTe, dtwopsi2F1112dnpar2;
0758     double twopsi2I1112,  dtwopsi2I1112dY, dtwopsi2I1112dTe, dtwopsi2I1112dnpar2;
0759     double dF72dY, dF72dTe, dF72dnpar2, dI72dY, dI72dTe, dI72dnpar2;
0760     double dF92dY, dF92dTe, dF92dnpar2, dI92dY, dI92dTe, dI92dnpar2;
0761     double F921, dF921dY, dF921dTe, dF921dnpar2, I921, dI921dY, dI921dTe, dI921dnpar2;
0762     double twopsi2F1122, dtwopsi2F1122dY, dtwopsi2F1122dTe, dtwopsi2F1122dnpar2;
0763     double twopsi2I1122, dtwopsi2I1122dY, dtwopsi2I1122dTe, dtwopsi2I1122dnpar2;
0764     double dF93dY, dF93dTe, dF93dnpar2, dI93dY, dI93dTe, dI93dnpar2;
0765     double delta_1, delta_1_1, delta_1_2, delta_2, delta_2_1, delta_2_2, delta_3, delta_3_1;
0766     double F5_1, dF5_1dz, muF5_1, dmuF5_1dY, dmuF5_1dTe, dmuF5_1dnpar2;
0767     double F7_1, dF7_1dz, dF7_1dY, dF7_1dTe, dF7_1dnpar2;
0768     double F9_1, dF9_1dz, dF9_1dY, dF9_1dTe, dF9_1dnpar2;
0769     double muF7_11, dmuF7_11dY, dmuF7_11dTe, dmuF7_11dnpar2;
0770     double F9_11, dF9_11dY, dF9_11dTe, dF9_11dnpar2;
0771     double twopsi2F9_12,  dtwopsi2F9_12dY,  dtwopsi2F9_12dTe,  dtwopsi2F9_12dnpar2;
0772     double twopsi2F11_12, dtwopsi2F11_12dY, dtwopsi2F11_12dTe, dtwopsi2F11_12dnpar2;
0773     double F7_2, dF7_2dz, dF7_2dY, dF7_2dTe, dF7_2dnpar2;
0774     double F9_2, dF9_2dz, dF9_2dY, dF9_2dTe, dF9_2dnpar2;
0775     double F9_21, dF9_21dY, dF9_21dTe, dF9_21dnpar2;
0776     double twopsi2F11_22, dtwopsi2F11_22dY, dtwopsi2F11_22dTe, dtwopsi2F11_22dnpar2;
0777     double F9_3, dF9_3dz, dF9_3dY, dF9_3dTe, dF9_3dnpar2;
0778     double coeff1, coeff2, coeff3, coeff4, coeff5, coeff6, coeff7, coeff8;
0779 
0780     double complex R0, dR0dX, dR0dY, dR0dTe, dR0dnpar2;
0781     double complex R1, dR1dX, dR1dY, dR1dTe, dR1dnpar2;
0782     double complex R2, dR2dX, dR2dY, dR2dTe, dR2dnpar2;
0783     double complex RR, dRdX , dRdY , dRdTe , dRdnpar2;
0784     double complex L0, dL0dX, dL0dY, dL0dTe, dL0dnpar2;
0785     double complex L1, dL1dX, dL1dY, dL1dTe, dL1dnpar2;
0786     double complex L2, dL2dX, dL2dY, dL2dTe, dL2dnpar2;
0787     double complex L , dLdX , dLdY , dLdTe , dLdnpar2;
0788     double complex P0, dP0dX, dP0dY, dP0dTe, dP0dnpar2;
0789     double complex P1, dP1dX, dP1dY, dP1dTe, dP1dnpar2;
0790     double complex P2, dP2dX, dP2dY, dP2dTe, dP2dnpar2;
0791     double complex P , dPdX , dPdY , dPdTe , dPdnpar2;
0792     double complex T1, dT1dX, dT1dY, dT1dTe, dT1dnpar2;
0793     double complex T2, dT2dX, dT2dY, dT2dTe, dT2dnpar2;
0794     double complex T,  dTdX,  dTdY,   dTdTe, dTdnpar2;
0795     double complex M0, dM0dX, dM0dY, dM0dTe, dM0dnpar2;
0796     double complex M1, dM1dX, dM1dY, dM1dTe, dM1dnpar2;
0797     double complex M, dMdX, dMdY, dMdTe, dMdnpar2;
0798     double complex N0, dN0dX, dN0dY, dN0dTe, dN0dnpar2;
0799     double complex N1, dN1dX, dN1dY, dN1dTe, dN1dnpar2;
0800     double complex N, dNdX,  dNdY,  dNdTe,  dNdnpar2;
0801     double complex S0, S1, S2, S;
0802     double complex Nperp2;
0803 
0804     
0805 
0806 /* Hot plasma, both modes FLR approximation up to (kperp*rhoe)^4
0807    Te is normalized Te => kTe/mc2 
0808    ne is normalized ne => X=e^2*ne/(me*epsilon0*omega0^2)
0809    B is normalized B => Y=e*B/(me*omega0)
0810     k is normalized k => k*clumn = refractive index 
0811 */  
0812 
0813     npar2 = Npar*Npar;
0814     te = betath[0]*betath[0];
0815     B = wcn[0];
0816     XX = wpn2[0];
0817     Y_1 = 1.0/B;
0818     Y_2 = Y_1*Y_1;
0819     mu = 1.0/te;
0820     halfmu = 0.5*mu;
0821     mu2 = mu*mu;
0822     delta1 = 1.0-B;
0823     delta2 = delta1-B;
0824     delta3 = delta2-B;
0825     halfnpar2 = 0.5*npar2;
0826     psi2 = mu*halfnpar2;
0827     twopsi2 = 2.0*psi2;
0828 
0829 /* Compute Shkarovsky functions together with its partial derivatives 
0830    for resonant parameters mu*(1.0-n*Y)  
0831    Non-resonant parameters (n<0) => asymptotic approximation 
0832 */
0833 
0834 
0835 /* Resonant terms N=0 */
0836     Shka(mu,psi2,&F50,&F70,&F90,&dF50dz,&dF70dz,&dF90dz,&dF30da,&dF50da,&dF70da,&dF90da,&dF110da);
0837 
0838     muF50 = mu*F50;
0839     dmuF50dTe = -mu2*(mu*(dF50dz+halfnpar2*dF50da)+F50);
0840     dmuF50dnpar2 = mu*halfmu*dF50da;
0841 
0842     dF70dTe = -mu2*(dF70dz+halfnpar2*dF70da);
0843     dF70dnpar2 = halfmu*dF70da; 
0844     dF90dTe = -mu2*(dF90dz+halfnpar2*dF90da);
0845     dF90dnpar2 = halfmu*dF90da; 
0846 
0847     muF701 = mu*dF70dz;
0848     dmuF701dTe = -mu2*(mu*(dF50da+halfnpar2*(dF70da-dF50da))+dF70dz);
0849     dmuF701dnpar2 = mu*halfmu*(dF70da-dF50da);
0850 
0851     F901 = dF90dz;
0852     dF901dTe = -mu2*(dF70da+halfnpar2*(dF90da-dF70da));
0853     dF901dnpar2 = halfmu*(dF90da-dF70da);
0854  
0855     twopsi2muF702 = mu*twopsi2*dF50da;
0856     dtwopsi2muF702dTe = -mu2*mu*twopsi2*(dF50da-dF30da+halfnpar2*(dF70da-2.0*dF50da+dF30da)) - 2.0*mu*twopsi2muF702;
0857     dtwopsi2muF702dnpar2 = mu*twopsi2*halfmu*(dF70da-2.0*dF50da+dF30da)+mu2*dF50da;
0858 
0859     twopsi2F902 = twopsi2*dF70da;
0860     dtwopsi2F902dTe = -mu2*twopsi2*(dF70da-dF50da+halfnpar2*(dF90da-2.0*dF70da+dF50da)) - mu*twopsi2F902;
0861     dtwopsi2F902dnpar2 = twopsi2*halfmu*(dF90da-2.0*dF70da+dF50da)+mu*dF70da;
0862      
0863     twopsi2F1102 = twopsi2*dF90da;
0864     dtwopsi2F1102dTe = -mu2*twopsi2*(dF90da-dF70da+halfnpar2*(dF110da-2.0*dF90da+dF70da)) - mu*twopsi2F1102;
0865     dtwopsi2F1102dnpar2=twopsi2*halfmu*(dF110da-2.0*dF90da+dF70da)+mu*dF90da;
0866 
0867     if (mu < psi2) {
0868         ImShka(mu,psi2,&I50,&I70,&I90,&dI50dz,&dI70dz,&dI90dz,&dI30da,&dI50da,&dI70da,&dI90da,&dI110da);
0869         muI50 = mu*I50;
0870         dmuI50dTe = -mu2*(mu*(dI50dz+halfnpar2*dI50da)+I50);
0871         dmuI50dnpar2 = mu*halfmu*dI50da;
0872         dI70dTe = -mu2*(dI70dz+halfnpar2*dI70da);
0873         dI70dnpar2 = halfmu*dI70da; 
0874         dI90dTe = -mu2*(dI90dz+halfnpar2*dI90da);
0875         dI90dnpar2 = halfmu*dI90da;
0876         muI701 = mu*dI70dz;
0877         dmuI701dTe = -mu2*(mu*(dI50da+halfnpar2*(dI70da-dI50da))+dI70dz);
0878         dmuI701dnpar2 = mu*halfmu*(dI70da-dI50da);
0879         I901 = dI90dz;
0880         dI901dTe = -mu2*(dI70da+halfnpar2*(dI90da-dI70da));
0881         dI901dnpar2 = halfmu*(dI90da-dI70da);
0882         twopsi2muI702 = mu*twopsi2*dI50da;
0883         dtwopsi2muI702dTe = -mu2*mu*twopsi2*(dI50da-dI30da+halfnpar2*(dI70da-2.0*dI50da+dI30da)) - 2.0*mu*twopsi2muI702;
0884         dtwopsi2muI702dnpar2 = mu*twopsi2*halfmu*(dI70da-2.0*dI50da+dI30da)+mu2*dI50da;
0885         twopsi2I902 = twopsi2*dI70da;
0886         dtwopsi2I902dTe = -mu2*twopsi2*(dI70da-dI50da+halfnpar2*(dI90da-2.0*dI70da+dI50da)) - mu*twopsi2I902;
0887         dtwopsi2I902dnpar2 = twopsi2*halfmu*(dI90da-2.0*dI70da+dI50da)+mu*dI70da;
0888         twopsi2I1102 = twopsi2*dI90da;
0889         dtwopsi2I1102dTe = -mu2*twopsi2*(dI90da-dI70da+halfnpar2*(dI110da-2.0*dI90da+dI70da)) - mu*twopsi2I1102;
0890         dtwopsi2I1102dnpar2 = twopsi2*halfmu*(dI110da-2.0*dI90da+dI70da)+mu*dI90da;
0891     } else {
0892         I50 = 0.0;I70 = 0.0; I90 = 0.0; dI50dz = 0.0; dI70dz = 0.0; dI90dz = 0.0;
0893         dI30da = 0.0; dI50da = 0.0; dI70da = 0.0; dI90da = 0.0; dI110da = 0.0;
0894         muI50 = 0.0; dmuI50dTe = 0.0; dmuI50dnpar2 = 0.0;
0895         dI70dTe = 0.0; dI70dnpar2 = 0.0;
0896         dI90dTe = 0.0; dI90dnpar2 = 0.0;
0897         muI701 = 0.0; dmuI701dTe = 0.0; dmuI701dnpar2 = 0.0;
0898         I901 = 0.0;   dI901dTe = 0.0;   dI901dnpar2 = 0.0;
0899         twopsi2muI702 = 0.0; dtwopsi2muI702dTe = 0.0; dtwopsi2muI702dnpar2 = 0.0;
0900         twopsi2I902  = 0.0; dtwopsi2I902dTe  = 0.0; dtwopsi2I902dnpar2  = 0.0;
0901         twopsi2I1102 = 0.0; dtwopsi2I1102dTe = 0.0; dtwopsi2I1102dnpar2 = 0.0;
0902     }
0903 
0904 /* Resonant terms N=+1 */
0905     Shka(mu*delta1,psi2,&F51,&F71,&F91,&dF51dz,&dF71dz,&dF91dz,&pipo1,&dF51da,&dF71da,&dF91da,&dF111da);         
0906 
0907     muF51 = mu*F51;
0908     dmuF51dY = -mu2*dF51dz;
0909     dmuF51dTe = -mu2*(mu*(delta1*dF51dz + halfnpar2*dF51da)+F51);
0910     dmuF51dnpar2 = mu*halfmu*dF51da;
0911 
0912     dF71dY = -mu*dF71dz;
0913     dF71dTe = -mu2*(delta1*dF71dz + halfnpar2*dF71da);
0914     dF71dnpar2 = halfmu*dF71da;
0915 
0916     dF91dY = -mu*dF91dz;
0917     dF91dTe = -mu2*(delta1*dF91dz + halfnpar2*dF91da);
0918     dF91dnpar2 = halfmu*dF91da;
0919 
0920     muF711 = mu*dF71dz;
0921     dmuF711dY = -mu2*dF51da;
0922     dmuF711dTe = -mu2*(mu*(delta1*dF51da+halfnpar2*(dF71da-dF51da))+dF71dz);
0923     dmuF711dnpar2 = mu*halfmu*(dF71da-dF51da);
0924 
0925     F911 = dF91dz;
0926     dF911dY = -mu*dF71da;
0927     dF911dTe = -mu2*(delta1*dF71da+halfnpar2*(dF91da-dF71da));
0928     dF911dnpar2 = halfmu*(dF91da-dF71da); 
0929      
0930     twopsi2F912 = twopsi2*dF71da;
0931     dtwopsi2F912dY = -twopsi2*mu*(dF71da-dF51da);
0932     dtwopsi2F912dTe = -mu2*twopsi2*(delta1*(dF71da-dF51da) + halfnpar2*(dF91da-2.0*dF71da+dF51da)) - mu*twopsi2F912;
0933     dtwopsi2F912dnpar2 = twopsi2*halfmu*(dF91da-2.0*dF71da+dF51da)+mu*dF71da;
0934      
0935     twopsi2F1112 = twopsi2*dF91da;
0936     dtwopsi2F1112dY = -twopsi2*mu*(dF91da-dF71da);
0937     dtwopsi2F1112dTe = -mu2*twopsi2*(delta1*(dF91da-dF71da) + halfnpar2*(dF111da-2.0*dF91da+dF71da)) - mu*twopsi2F1112;
0938     dtwopsi2F1112dnpar2 = twopsi2*halfmu*(dF111da-2.0*dF91da+dF71da)+mu*dF91da;
0939 
0940     if (mu*delta1else {
0969         I51 = 0.0;
0970     I71 = 0.0;
0971     I91 = 0.0;
0972     dI51dz = 0.0;
0973     dI71dz = 0.0;
0974     dI91dz = 0.0;
0975         dI51da = 0.0;
0976     dI71da = 0.0;
0977     dI91da = 0.0;
0978     dI111da = 0.0;
0979         muI51 = 0.0;
0980     dmuI51dY = 0.0;
0981     dmuI51dTe = 0.0;
0982     dmuI51dnpar2 = 0.0;
0983         dI71dY = 0.0;
0984     dI71dTe = 0.0;
0985     dI71dnpar2 = 0.0;
0986         dI91dY = 0.0;
0987     dI91dTe = 0.0;
0988     dI91dnpar2 = 0.0;
0989         muI711 = 0.0;
0990     dmuI711dY = 0.0;
0991     dmuI711dTe = 0.0;
0992     dmuI711dnpar2 = 0.0;
0993         I911 = 0.0;
0994     dI911dY = 0.0;
0995     dI911dTe = 0.0;
0996     dI911dnpar2 = 0.0;
0997         twopsi2I912  = 0.0;
0998     dtwopsi2I912dY  = 0.0;
0999     dtwopsi2I912dTe  = 0.0;
1000     dtwopsi2I912dnpar2  = 0.0;
1001         twopsi2I1112 = 0.0;
1002     dtwopsi2I1112dY = 0.0;
1003     dtwopsi2I1112dTe = 0.0;
1004     dtwopsi2I1112dnpar2 = 0.0;
1005     }
1006 
1007 /* Resonant terms N = +2 */
1008     Shka(mu*delta2,psi2,&pipo1,&F72,&F92,&pipo1,&dF72dz,&dF92dz,&pipo1,&pipo1,&dF72da,&dF92da,&dF112da);         
1009 
1010     dF72dY = -2.0*mu*dF72dz;
1011     dF72dTe = -mu2*(delta2*dF72dz+halfnpar2*dF72da);
1012     dF72dnpar2 = halfmu*dF72da;
1013 
1014     dF92dY = -2.0*mu*dF92dz;
1015     dF92dTe = -mu2*(delta2*dF92dz+halfnpar2*dF92da);
1016     dF92dnpar2 = halfmu*dF92da;
1017 
1018     F921 = dF92dz;
1019     dF921dY = -2.0*mu*dF72da;
1020     dF921dTe = -mu2*(delta2*dF72da+halfnpar2*(dF92da-dF72da));
1021     dF921dnpar2 = halfmu*(dF92da-dF72da);
1022 
1023     twopsi2F1122 = twopsi2*dF92da;
1024     dtwopsi2F1122dY = -2.0*twopsi2*mu*(dF92da-dF72da);
1025     dtwopsi2F1122dTe = -mu2*twopsi2*(delta2*(dF92da-dF72da) + halfnpar2*(dF112da-2.0*dF92da+dF72da)) - mu*twopsi2F1122;
1026     dtwopsi2F1122dnpar2 = twopsi2*halfmu*(dF112da-2.0*dF92da+dF72da)+mu*dF92da; 
1027 
1028     if (mu*delta2else {
1046         I72 = 0.0;
1047     I92 = 0.0;
1048     dI72dz = 0.0;
1049     dI92dz = 0.0;
1050         dI72da = 0.0;
1051     dI92da = 0.0;
1052     dI112da = 0.0;         
1053         dI72dY = 0.0;
1054     dI72dTe = 0.0;
1055     dI72dnpar2 = 0.0;
1056         dI92dY = 0.0;
1057     dI92dTe = 0.0;
1058     dI92dnpar2 = 0.0;
1059         I921 = 0.0;
1060     dI921dY = 0.0;
1061     dI921dTe = 0.0;
1062     dI921dnpar2 = 0.0;
1063         twopsi2I1122 = 0.0;
1064     dtwopsi2I1122dY = 0.0;
1065     dtwopsi2I1122dTe = 0.0;
1066     dtwopsi2I1122dnpar2 = 0.0;
1067     }
1068 
1069 /* Resonant terms N = +3 */
1070     
1071     Shka(mu*delta3,psi2,&pipo1,&pipo1,&F93,&pipo1,&pipo1,&dF93dz,&pipo1,&pipo1,&pipo1,&dF93da,&pipo1);         
1072 
1073     dF93dY = -3.0*mu*dF93dz;
1074     dF93dTe = -mu2*(delta3*dF93dz+halfnpar2*dF93da);
1075     dF93dnpar2 = halfmu*dF93da;
1076 
1077     if (mu*delta3else {
1084         I93 = 0.0; dI93dz = 0.0;
1085         dI93da = 0.0;
1086     dI93dY = 0.0;
1087     dI93dTe = 0.0;
1088     dI93dnpar2 = 0.0;
1089     }
1090 
1091          
1092 /* Non-resonant parameters (n<0)  = > asymptotic approximation for Shkarovski functions */
1093 
1094     delta_1 = 1.0+B;
1095     delta_2 = delta_1+B;
1096     delta_3 = delta_2+B;
1097     delta_1_1 = 1.0/delta_1;
1098     delta_1_2 = delta_1_1*delta_1_1;
1099     delta_2_1 = 1.0/delta_2;
1100     delta_2_2 = delta_2_1*delta_2_1;
1101     delta_3_1 = 1.0/delta_3;
1102 
1103     muF5_1 = delta_1_1;
1104     dmuF5_1dY = -delta_1_2;
1105     dmuF5_1dTe = 0.0;
1106     dmuF5_1dnpar2 = 0.0;
1107 
1108     F7_1 = te*delta_1_1;
1109     dF7_1dY = -te*delta_1_2;
1110     dF7_1dTe = delta_1_1;
1111     dF7_1dnpar2 = 0.0;
1112 
1113     F9_1 = F7_1; dF9_1dY = dF7_1dY; dF9_1dTe = dF7_1dTe; dF9_1dnpar2 = 0.0;
1114 
1115     muF7_11 = -te*delta_1_2;
1116     dmuF7_11dY = -2.0*muF7_11*delta_1_1;
1117     dmuF7_11dTe = mu*muF7_11;
1118     dmuF7_11dnpar2 = 0.0;
1119 
1120     F9_11 = -te*te*delta_1_2;
1121     dF9_11dY = -2.0*F9_11*delta_1_1;
1122     dF9_11dTe = 2.0*mu*F9_11;
1123     dF9_11dnpar2 = 0.0; 
1124 
1125     dtwopsi2F9_12dnpar2 = 2.0*te*te*delta_1_2*delta_1_1;
1126     twopsi2F9_12 = npar2*dtwopsi2F9_12dnpar2;
1127     dtwopsi2F9_12dY = -3.0*delta_1_1*twopsi2F9_12;
1128     dtwopsi2F9_12dTe = 2.0*mu*twopsi2F9_12;
1129 
1130     twopsi2F11_12 = twopsi2F9_12; dtwopsi2F11_12dY = dtwopsi2F9_12dY; 
1131     dtwopsi2F11_12dTe = dtwopsi2F9_12dTe; dtwopsi2F11_12dnpar2 = dtwopsi2F9_12dnpar2; 
1132      
1133     F7_2 = te*delta_2_1;
1134     dF7_2dY = -2.0*te*delta_2_2;
1135     dF7_2dTe = delta_2_1;
1136     dF7_2dnpar2 = 0.0;
1137 
1138     F9_2 = F7_2; dF9_2dY = dF7_2dY; dF9_2dTe = dF7_2dTe; dF9_2dnpar2 = 0.0;
1139 
1140     F9_21 = -te*te*delta_2_2;
1141     dF9_21dY = -4.0*F9_21*delta_2_1;
1142     dF9_21dTe = 2.0*mu*F9_21;
1143     dF9_21dnpar2 = 0.0;
1144 
1145     dtwopsi2F11_22dnpar2 = 2.0*te*te*delta_2_2*delta_2_1;
1146     twopsi2F11_22 = npar2*dtwopsi2F11_22dnpar2;
1147     dtwopsi2F11_22dY = -6.0*delta_2_1*twopsi2F11_22;
1148     dtwopsi2F11_22dTe = 2.0*mu*twopsi2F11_22; 
1149 
1150     F9_3 = te*delta_3_1;
1151     dF9_3dY = -3.0*F9_3*delta_3_1;
1152     dF9_3dTe = delta_3_1;
1153     dF9_3dnpar2 = 0.0;
1154 
1155 /* Compute the weakly relativistic tensor elements in rotating frame.
1156    FLR approxomation up to (kperp*rhoe)^4, 
1157    Tensor is developped up to the +/- second harmonic of Wce.
1158    I.P. Shkarovsky, J. Plasma Phys. (1986) vol. 35 part2, pp319-331.
1159     = > valid for O1, X1, O2, X2, X3 ECRH scenarios */
1160 
1161     coeff1 = 0.5*Y_2;
1162     coeff2 = XX*coeff1;
1163     coeff4 = XX*Y_1;
1164     coeff3 = mu*coeff4;
1165     coeff5 = 0.5*te*coeff1*coeff1; /*  = te/(8*Y^4) */
1166     coeff6 = XX*coeff5;
1167     coeff7 = 0.25*Y_2*Y_1;
1168     coeff8 = XX*coeff7;
1169 
1170     dR0dX = -muF51-I*muI51;
1171     R0 = XX*dR0dX;
1172     dR0dY = -XX*(dmuF51dY+I*dmuI51dY);
1173     dR0dTe =  -XX*(dmuF51dTe+I*dmuI51dTe);
1174     dR0dnpar2 = -XX*(dmuF51dnpar2+I*dmuI51dnpar2);/* ligne a tester YP */
1175 
1176     dR1dX = coeff1*((3.0*F71-F7_1-2.0*F72)+I*(3.0*I71 - 2.0*I72));
1177     R1 = XX*dR1dX;
1178     dR1dY = -2.0*R1*Y_1+coeff2*((3.0*dF71dY-dF7_1dY-2.0*dF72dY) + I*(3.0*dI71dY-2.0*dI72dY));
1179     dR1dTe = coeff2*((3.0*dF71dTe - dF7_1dTe - 2.0*dF72dTe) + I*(3.0*dI71dTe - 2.0*dI72dTe));
1180     dR1dnpar2 = coeff2*((3.0*dF71dTe-dF7_1dTe-2.0*dF72dTe)+I*(3.0*dI71dTe-2.0*dI72dTe));
1181     dR2dX = coeff5*((5.0*F9_1-10.0*F91-2.0*F9_2-3.0*F93) + I*(-10.0*I91-3.0*I93));
1182     R2 = XX*dR2dX;
1183     dR2dY = -4.0*R2*Y_1+coeff6*((5.0*dF9_1dY-10.0*dF91dY-2.0*dF9_2dY-3.0*dF93dY) + I*(-10.0*dI91dY-3.0*dI93dY));
1184     dR2dTe = mu*R2+coeff6*((5.0*dF9_1dTe-10.0*dF91dTe-2.0*dF9_2dTe-3.0*dF93dTe) + I*(-10.0*dI91dTe-3.0*dI93dTe));
1185     dR2dnpar2 = coeff6*((5.0*dF9_1dnpar2-10.0*dF91dnpar2-2.0*dF9_2dnpar2-3.0*dF93dnpar2) + I*(-10.0*dI91dnpar2-3.0*dI93dnpar2));/*  sign Im term to be checked YP ******/
1186 
1187     dL0dX = -muF5_1;
1188     L0 = -XX*muF5_1;
1189     dL0dY = -XX*dmuF5_1dY;
1190     dL0dTe = -XX*dmuF5_1dTe;
1191     dL0dnpar2 = -XX*dmuF5_1dnpar2;
1192 
1193     dL1dX = coeff1*((3.0*F7_1-F71-2.0*F7_2) + I*(-I71));
1194     L1 = XX*dL1dX;
1195     dL1dY = -2.0*L1*Y_1+coeff2*((3.0*dF7_1dY-dF71dY-2.0*dF7_2dY) + I*(-dI71dY));
1196     dL1dTe = coeff2*((3.0*dF7_1dTe-dF71dTe-2.0*dF7_2dTe) +I*(-dI71dTe));
1197     dL1dnpar2 = coeff2*((3.0*dF7_1dnpar2-dF71dnpar2-2.0*dF7_2dnpar2) + I*(-dI71dnpar2)); 
1198 
1199     dL2dX = coeff5*((5.0*F91-10.0*F9_1-2.0*F92-3.0*F9_3) + I*(5.0*I91-2.0*I92));
1200     L2 = XX*dL2dX;
1201     dL2dY = -4.0*L2*Y_1+coeff6*((5.0*dF91dY-10.0*dF9_1dY-2.0*dF92dY-3.0*dF9_3dY) + I*(5.0*dI91dY-2.0*dI92dY));
1202     dL2dTe = mu*L2+coeff6*((5.0*dF91dTe-10.0*dF9_1dTe-2.0*dF92dTe-3.0*dF9_3dTe) + I*(5.0*dI91dTe-2.0*dI92dTe));
1203     dL2dnpar2 = coeff6*((5.0*dF91dnpar2-10.0*dF9_1dnpar2 -2.0*dF92dnpar2-3.0*dF9_3dnpar2) + I*(5.0*dI91dnpar2-2.0*dI92dnpar2));
1204 
1205     dP0dX = (-muF50-twopsi2muF702)+I*(-muI50-twopsi2muI702);
1206     P0 = XX*dP0dX;
1207     dP0dY = 0.0; 
1208     dP0dTe = -XX*((dmuF50dTe+dtwopsi2muF702dTe) + I*(dmuI50dTe+dtwopsi2muI702dTe));
1209     dP0dnpar2 = -XX*((dmuF50dnpar2 +dtwopsi2muF702dnpar2) + I*(dmuI50dnpar2 +dtwopsi2muI702dnpar2));
1210 
1211     dP1dX = coeff1*((2.0*(F70+twopsi2F902)-(F71+twopsi2F912)-(F7_1+twopsi2F9_12)) + I*(2.0*(I70+twopsi2I902)-(I71+twopsi2I912)));                            
1212     P1 = XX*dP1dX;
1213     dP1dY = -2.0*coeff4*dP1dX-coeff2*((dF71dY+dtwopsi2F912dY+dF7_1dY+dtwopsi2F9_12dY) + I*(dI71dY+dtwopsi2I912dY));
1214     dP1dTe = coeff2*((2.0*(dF70dTe+dtwopsi2F902dTe)-(dF71dTe+dtwopsi2F912dTe)-(dF7_1dTe+dtwopsi2F9_12dTe)) + I*(2.0*(dI70dTe+dtwopsi2I902dTe)-(dI71dTe+dtwopsi2I912dTe))); 
1215     dP1dnpar2 = coeff2*((2.0*(dF70dnpar2+dtwopsi2F902dnpar2)-(dF71dnpar2+dtwopsi2F912dnpar2)-(dF7_1dnpar2+dtwopsi2F9_12dnpar2)) + I*(2.0*(dI70dnpar2+dtwopsi2I902dnpar2)-(dI71dnpar2+dtwopsi2I912dnpar2)                                  ));
1216 
1217     dP2dX = coeff5*((-6.0*(F90+twopsi2F1102)+4.0*(F91+twopsi2F1112+F9_1+twopsi2F11_12)-(F92+twopsi2F1122+F9_2+twopsi2F11_22)) +I*(-6.0*(I90+twopsi2I1102)+4.0*(I91+twopsi2I1112)-(I92+twopsi2I1122                   )));
1218     P2 = XX*dP2dX;
1219     dP2dY = -4.0*P2*Y_1+coeff6*((+4.0*(dF91dY+dtwopsi2F1112dY+dF9_1dY+dtwopsi2F11_12dY) - (dF92dY+dtwopsi2F1122dY+dF9_2dY+dtwopsi2F11_22dY)) + I*(+4.0*(dI91dY+dtwopsi2I1112dY) - (dI92dY+dtwopsi2I1122dY                         )));
1220     dP2dTe = mu*P2+coeff6*((-6.0*(dF90dTe+dtwopsi2F1102dTe)+4.0*(dF91dTe+dtwopsi2F1112dTe+dF9_1dTe+dtwopsi2F11_12dTe)-(dF92dTe+dtwopsi2F1122dTe+dF9_2dTe+dtwopsi2F11_22dTe))+I*(-6.0*(dI90dTe+dtwopsi2I1102dTe)+4.0*(dI91dTe+dtwopsi2I1112dTe)-(dI92dTe+dtwopsi2I1122dTe)));
1221     dP2dnpar2 = coeff6*((-6.0*(dF90dnpar2+dtwopsi2F1102dnpar2)+4.0*(dF91dnpar2+dtwopsi2F1112dnpar2+dF9_1dnpar2+dtwopsi2F11_12dnpar2)-(dF92dnpar2+dtwopsi2F1122dnpar2+dF9_2dnpar2+dtwopsi2F11_22dnpar2))+I*(-6.0*(dI90dnpar2+dtwopsi2I1102dnpar2)+4.0*(dI91dnpar2+dtwopsi2I1112dnpar2)-(dI92dnpar2+dtwopsi2I1122dnpar2)));
1222 
1223     dT1dX = coeff1*((2.0*F70-F71-F7_1)+I*(2.0*I70-I71));
1224     T1 = XX*dT1dX;
1225     dT1dY = -2.0*coeff4*dT1dX-coeff2*((dF71dY+dF7_1dY) + I*(dI71dY));/*  sign Im term to be checked YP ******/
1226     dT1dTe = coeff2*((2.0*dF70dTe-dF71dTe-dF7_1dTe)+I*(2.0*dI70dTe-dI71dTe));
1227     dT1dnpar2 = coeff2*((2.0*dF70dnpar2 -dF71dnpar2-dF7_1dnpar2)+I*(2.0*dI70dnpar2-dI71dnpar2));
1228 
1229     dT2dX = coeff5*((16.0*(F91+F9_1)+4.0*(F92+F9_2)-12.0*F90)+I*(16.0*(I91)+4.0*(I92)-12.0*I90));
1230     T2 = XX*dT2dX;
1231     dT2dY = -4.0*Y_1*T2+coeff6*((16.0*(dF91dY+dF9_1dY)+4.0*(dF92dY+dF9_2dY)) +I*(16.0* dI91dY+4.0* dI92dY));
1232     dT2dTe = mu*T2+coeff6*((16.0*(dF91dTe+dF9_1dTe)+4.0*(dF92dTe+dF9_2dTe)-12.0*dF90dTe) +I*(16.0* dI91dTe+4.0* dI92dTe) -12.0*dI90dTe);
1233     dT2dnpar2 = coeff6*((16.0*(dF91dnpar2+dF9_1dnpar2)+4.0*(dF92dnpar2+dF9_2dnpar2)-12.0*dF90dnpar2)+I*(16.0* dI91dnpar2+4.0* dI92dnpar2-12.0*dI90dnpar2));
1234 
1235 
1236 /* L. Colas 11-02-2004 */
1237 /* the sign of dM0dX is not clearly established */
1238          
1239     dM0dX = Y_1*(-muF7_11+muF701+I*muI701);
1240     M0 = XX*dM0dX;
1241     dM0dY = coeff4*(-dmuF7_11dY-dM0dX);
1242     dM0dTe = coeff4*(-dmuF7_11dTe+dmuF701dTe+I*dmuI701dTe);
1243     dM0dnpar2 = coeff4*(-dmuF7_11dnpar2+dmuF701dnpar2+I*dmuI701dnpar2);
1244 
1245 
1246 /* The sign of M1 is not clear */
1247 
1248     dM1dX = coeff7*((2.0*F911+6.0*F9_11+2.0*F9_21-3.0*F901)+I*(2.0*I911-3.0*I901));
1249     M1 = XX*dM1dX;
1250     dM1dY = -3.0*M1*Y_1+coeff8*((2.0*dF911dY+6.0*dF9_11dY+2.0*dF9_21dY)+I*(2.0*dI911dY));
1251     dM1dTe = coeff8*((2.0*dF911dTe+6.0*dF9_11dTe+2.0*dF9_21dTe   -3.0*dF901dTe)+I*(2.0*dI911dTe-3.0*dI901dTe));
1252     dM1dnpar2 = coeff8*((2.0*dF911dnpar2+6.0*dF9_11dnpar2+2.0*dF9_21dnpar2-3.0*dF901dnpar2)+I*(2.0*dI911dnpar2-3.0*dI901dnpar2)); 
1253 
1254 
1255 /* L. Colas 11-02-2004 */
1256 /* the sign of dN0dX is not clearly established */
1257 
1258     dN0dX = Y_1*((muF711-muF701)+I*(muI711-muI701));
1259     N0 = XX*dN0dX;
1260     dN0dY = coeff4*(dmuF711dY-dN0dX+I*dmuI711dY);
1261     dN0dTe = coeff4*((dmuF711dTe-dmuF701dTe)+I*(dmuI711dTe-dmuI701dTe));
1262     dN0dnpar2 = coeff4*(dmuF711dnpar2-dmuF701dnpar2) + I*(dmuI711dnpar2-dmuI701dnpar2);/*  ligne � tester YP */
1263 
1264 
1265 /* The sign of N1 is not clear */
1266 
1267     dN1dX = coeff7*((-2.0*F9_11-6.0*F911-2.0*F921+3.0*F901)+I*(-6.0*I911-2.0*I921+3.0*I901));
1268     N1 = XX*dN1dX;
1269     dN1dY = -3.0*N1*Y_1+coeff8*((-2.0*dF9_11dY-6.0*dF911dY-2.0*dF921dY) +I*(-6.0*dI911dY-2.0*dI921dY));
1270     dN1dTe = coeff8*((-2.0*dF9_11dTe-6.0*dF911dTe-2.0*dF921dTe+3.0*dF901dTe) +I*(-6.0*dI911dTe-2.0*dI921dTe+3.0*dI901dTe));
1271     dN1dnpar2 = coeff8*((-2.0*dF9_11dnpar2-6.0*dF911dnpar2-2.0*dF921dnpar2+3.0*dF901dnpar2)+I*(-6.0*dI911dnpar2-2.0*dI921dnpar2+3.0*dI901dnpar2)); 
1272 
1273          
1274     Nperp2 = Nperp*Nperp;
1275     RR = R0+Nperp2*(R1+Nperp2*R2);
1276     L = L0+Nperp2*(L1+Nperp2*L2);
1277     S = 0.5*(RR+L);
1278     P = P0+Nperp2*(P1+Nperp2*P2);
1279     T = Nperp2*(T1+Nperp2*T2);
1280     M = M0+Nperp2*M1;
1281     N = N0+Nperp2*N1;
1282    
1283     dRdX = dR0dX+Nperp2*(dR1dX+Nperp2*dR2dX);
1284     dLdX = dL0dX+Nperp2*(dL1dX+Nperp2*dL2dX);
1285     dPdX = dP0dX+Nperp2*(dP1dX+Nperp2*dP2dX);
1286     dTdX =       Nperp2*(dT1dX+Nperp2*dT2dX);
1287     dMdX = dM0dX+Nperp2*dM1dX;
1288     dNdX = dN0dX+Nperp2*dN1dX;
1289 
1290     dRdY = dR0dY+Nperp2*(dR1dY+Nperp2*dR2dY);
1291     dLdY = dL0dY+Nperp2*(dL1dY+Nperp2*dL2dY);
1292     dPdY = dP0dY+Nperp2*(dP1dY+Nperp2*dP2dY);
1293     dTdY =       Nperp2*(dT1dY+Nperp2*dT2dY);
1294     dMdY = dM0dY+Nperp2*dM1dY;
1295     dNdY = dN0dY+Nperp2*dN1dY;
1296         
1297     dRdTe = dR0dTe+Nperp2*(dR1dTe+Nperp2*dR2dTe);
1298     dLdTe = dL0dTe+Nperp2*(dL1dTe+Nperp2*dL2dTe);
1299     dPdTe = dP0dTe+Nperp2*(dP1dTe+Nperp2*dP2dTe);
1300     dTdTe =        Nperp2*(dT1dTe+Nperp2*dT2dTe);
1301     dMdTe = dM0dTe+Nperp2*dM1dTe;
1302     dNdTe = dN0dTe+Nperp2*dN1dTe;
1303 
1304     dRdnpar2 = dR0dnpar2+Nperp2*(dR1dnpar2+Nperp2*dR2dnpar2);
1305     dLdnpar2 = dL0dnpar2+Nperp2*(dL1dnpar2+Nperp2*dL2dnpar2);
1306     dPdnpar2 = dP0dnpar2+Nperp2*(dP1dnpar2+Nperp2*dP2dnpar2);
1307     dTdnpar2 =           Nperp2*(dT1dnpar2+Nperp2*dT2dnpar2);
1308     dMdnpar2 = dM0dnpar2+Nperp2*dM1dnpar2;
1309     dNdnpar2 = dN0dnpar2+Nperp2*dN1dnpar2;
1310 
1311 
1312 /* Dielectric tensor for electrons, STIX frame */
1313 /* Marco BRAMBILLA, "Kinetic theory of plasma waves" 
1314     equation (24.1) */
1315    
1316     Xxx[0] = S;
1317     Xxy[0] = 0.5*I*(L-RR);
1318     Xxz[0] = 0.5*Npar*Nperp*(M+N);
1319     Xyy[0] = S-2.0*T;
1320     Xyz[0] = 0.5*I*Npar*Nperp*(M-N);
1321     Xzz[0] = P;
1322 
1323 /* Ionic susceptibility is negligible at high frequencies */
1324         
1325     for (is = 1;is < ns + 1;is++) {
1326         Xxx[is] = 0.0;
1327         Xxy[is] = 0.0;
1328         Xxz[is] = 0.0;
1329         Xyy[is] = 0.0;
1330         Xyz[is] = 0.0;
1331         Xzz[is] = 0.0;
1332     }
1333  
1334 /* derivatives of the dielectric tensor for electrons, STIX frame */
1335 
1336     dXxxdNperp[0] = Nperp*(R1+L1+2.0*Nperp2*(R2+L2));
1337     dXxydNperp[0] = I*Nperp*(L1-R1+2.0*Nperp2*(L2-R2));
1338     dXxzdNperp[0] = Npar*(0.5*(M+N)+Nperp2*(M1+N1));
1339     dXyydNperp[0] = dXxxdNperp[0]-4.0*Nperp*(T1+2.0*Nperp2*T2);
1340     dXyzdNperp[0] = I*Npar*(0.5*(M-N)+Nperp2*(M1-N1));
1341     dXzzdNperp[0] = 2.0*Nperp*(P1+2.0*Nperp2*P2);        
1342 
1343     dXxxdNpar[0] = Npar*(dRdnpar2+dLdnpar2);
1344     dXxydNpar[0] = Npar*I*(dLdnpar2-dRdnpar2);
1345     dXxzdNpar[0] = npar2*Nperp*(dMdnpar2+dNdnpar2)+0.5*Nperp*(M+N);
1346     dXyydNpar[0] = dXxxdNpar[0]-4.0*Npar*dTdnpar2;
1347     dXyzdNpar[0] = I*npar2*Nperp*(dMdnpar2-dNdnpar2)+0.5*I*Nperp*(M-N);;
1348     dXzzdNpar[0] = 2.0*Npar*dPdnpar2;
1349 
1350     dXxxdbetath[0] = betath[0]*(dRdTe+dLdTe);
1351     dXxydbetath[0] = betath[0]*I*(dLdTe-dRdTe);
1352     dXxzdbetath[0] = betath[0]*Npar*Nperp*(dMdTe+dNdTe);
1353     dXyydbetath[0] = 2.0*betath[0]*(0.5*(dRdTe+dLdTe)-2.0*dTdTe);
1354     dXyzdbetath[0] = betath[0]*I*Npar*Nperp*(dMdTe-dNdTe);
1355     dXzzdbetath[0] = 2.0*betath[0]*dPdTe;
1356 
1357     dXxxdwcn[0] = 0.5*(dRdY+dLdY);
1358     dXxydwcn[0] = 0.5*I*(dLdY-dRdY);
1359     dXxzdwcn[0] = 0.5*Npar*Nperp*(dMdY+dNdY);
1360     dXyydwcn[0] = 0.5*(dRdY+dLdY)-2.0*dTdY;
1361     dXyzdwcn[0] = 0.5*I*Npar*Nperp*(dMdY-dNdY);
1362     dXzzdwcn[0] = dPdY;
1363  
1364     dXxxdwpn2[0] = 0.5*(dRdX+dLdX);
1365     dXxydwpn2[0] = 0.5*I*(dLdX-dRdX);
1366     dXxzdwpn2[0] = 0.5*Npar*Nperp*(dMdX+dNdX);
1367     dXyydwpn2[0] = 0.5*(dRdX+dLdX)-2.0*dTdX;
1368     dXyzdwpn2[0] = 0.5*I*Npar*Nperp*(dMdX-dNdX);
1369     dXzzdwpn2[0] = dPdX;    
1370 
1371 /* Ionic susceptibility is negligible at high frequencies */
1372         
1373     for (is = 1;is < ns + 1; is++) {
1374         dXxxdNperp[is] = 0.0;
1375         dXxydNperp[is] = 0.0;
1376         dXxzdNperp[is] = 0.0;
1377         dXyydNperp[is] = 0.0;
1378         dXyzdNperp[is] = 0.0;
1379         dXzzdNperp[is] = 0.0;
1380         dXxxdNpar[is] = 0.0;
1381         dXxydNpar[is] = 0.0;
1382         dXxzdNpar[is] = 0.0;
1383         dXyydNpar[is] = 0.0;
1384         dXyzdNpar[is] = 0.0;
1385         dXzzdNpar[is] = 0.0;
1386         dXxxdbetath[is] = 0.0;
1387         dXxydbetath[is] = 0.0;
1388         dXxzdbetath[is] = 0.0;
1389         dXyydbetath[is] = 0.0;
1390         dXyzdbetath[is] = 0.0;
1391         dXzzdbetath[is] = 0.0;
1392         dXxxdwcn[is] = 0.0;
1393         dXxydwcn[is] = 0.0;
1394         dXxzdwcn[is] = 0.0;
1395         dXyydwcn[is] = 0.0;
1396         dXyzdwcn[is] = 0.0;
1397         dXzzdwcn[is] = 0.0;    
1398         dXxxdwpn2[is] = 0.0;
1399         dXxydwpn2[is] = 0.0;
1400         dXxzdwpn2[is] = 0.0;
1401         dXyydwpn2[is] =0.0;
1402         dXyzdwpn2[is] = 0.0;
1403         dXzzdwpn2[is] = 0.0;    
1404     }
1405     buildXstruc(X,dNperpdx,Xxx,Xxy,Xxz,Xyy,Xyz,Xzz,dXxxdNperp,dXxydNperp,dXxzdNperp,dXyydNperp,dXyzdNperp,dXzzdNperp,dXxxdNpar,dXxydNpar,dXxzdNpar,dXyydNpar,dXyzdNpar,dXzzdNpar,dXxxdbetath,dXxydbetath,dXxzdbetath,dXyydbetath,dXyzdbetath,dXzzdbetath,dXxxdwcn,dXxydwcn,dXxzdwcn,dXyydwcn,dXyzdwcn,dXzzdwcn,dXxxdwpn2,dXxydwpn2,dXxzdwpn2,dXyydwpn2,dXyzdwpn2,dXzzdwpn2,ns,dNpardx,dwcndx,dwpn2dx,dbetathdx);
1406 }
1407 
1408 
1409 
1410 void Shka(z,a,RF5,RF7,RF9,dRF5dz,dRF7dz,dRF9dz,dRF3da,dRF5da,dRF7da,dRF9da,dRF11da)
1411 double z,a;
1412 double *RF5,*RF7,*RF9,*dRF5dz,*dRF7dz,*dRF9dz,*dRF3da,*dRF5da,*dRF7da,*dRF9da,*dRF11da;
1413 {
1414 
1415 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1416 %% Real part of Shkarovsky functions F_q(z,a)
1417 %% for q = 1.5 2.5 3.5, along with their partial derivatives
1418 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1419 %% P.A. Robinson, J. Math. Phys. vol. 27, N�5, May 1986
1420 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1421 */
1422 
1423     double S0, S1, S2, S3, S4, S5, S6;
1424     double dS1dz, dS2dz, dS3dz, dS4dz, dS5dz, dS6dz; 
1425     double dS1da, dS2da, dS3da, dS4da, dS5da;
1426     double F1, F2, F3, F4, F5, F6, dF1, dF2, dF3, dF4, dF5, dF6;
1427     double d2F1, d2F2, d2F3, d2F4, d2F5;
1428     double qmax, limit, expo, oldcoeff, az1, terme, terme1, nmax, z_a, coeff;
1429     double phi2, phi, psi, rZ1, rZ2;
1430     double Cqj[1000];
1431     double daw;
1432     
1433     int nb, n, compteur;
1434     static double sqrttwoe_6 = 1.414e-6;
1435 /* LC 08/10/2006 : Zfried replaced with ZfriedConte */
1436     double complex Z1, Z2;
1437 
1438 
1439 /* %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1440 %% Asymptotic form
1441 %% P.A. Robinson, J. Math. Phys. vol. 27, N�5, May 1986
1442 %% formula (35) & (36)
1443 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% */
1444 
1445     qmax = 6.0;
1446     limit = 10*sqrt(a*a+qmax*qmax);
1447     nb = (int) limit;
1448     nb = ((nb) > (1000) ? (1000) : (nb)); /* nb is limited to 1000 */
1449 
1450     if ((fabs(z-a)>10.0*limit)&(a>sqrttwoe_6)) {
1451         for (n = 0;n/* !! = j+1 in Robinson formula !! */
1456             az1 = 1.0/(a-z);
1457             terme = az1;
1458         /* S1 = 0.0; */
1459         S2 = 0.0;
1460         S3 = 0.0; 
1461         S4 = 0.0; /* S5 = 0.0; S6 = 0.0; */
1462         /* dS1dz = 0.0 ; */
1463         dS2dz = 0.0;
1464         dS3dz = 0.0;
1465         dS4dz = 0.0; /* dS5dz = 0.0; dS6dz = 0.0; */
1466         dS1da = 0.0;
1467         dS2da = 0.0;
1468         dS3da = 0.0; 
1469         dS4da = 0.0;
1470         dS5da = 0.0;
1471       
1472         while (((Cqj[0]/oldcoeff)0.0)) {
1473 
1474 /* S1- = (Cqj[0]*terme); */
1475                 S2 -= (Cqj[1]*terme);
1476                 S3 -= (Cqj[2]*terme);
1477                 S4 -= (Cqj[3]*terme);
1478 /* S5- = (Cqj[4]*terme); */        
1479 /* S6- = (Cqj[5]*terme); */
1480                 terme *= az1;
1481                 terme1 = -compteur*terme;      
1482 /* dS1dz+ = (Cqj[0]*terme1); */
1483                 dS2dz += (Cqj[1]*terme1);
1484                 dS3dz += (Cqj[2]*terme1);
1485                 dS4dz += (Cqj[3]*terme1);
1486 /* dS5dz+ = (Cqj[4]*terme1);*/
1487 /* dS6dz+ = (Cqj[5]*terme1);*/
1488                 terme1 *= ((compteur+1.0)*az1);        
1489                 dS1da += (Cqj[1]*terme1);
1490                 dS2da += (Cqj[2]*terme1);
1491                 dS3da += (Cqj[3]*terme1);
1492                 dS4da += (Cqj[4]*terme1);
1493                 dS5da += (Cqj[5]*terme1);
1494                 oldcoeff = Cqj[0];
1495     /* %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1496     %% Update coefficients
1497     %% P.A. Robinson, J. Math. Phys. vol. 27, N�5, May 1986
1498     %% formula (36)
1499     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% */
1500             for (n = 0;n/* end while */
1505    
1506     }   /* endif */
1507 
1508     /* %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1509     %% value for a = 0
1510     %% P.A. Robinson, J. Math. Phys. vol. 27, N�5, May 1986
1511     %% formula (3)
1512     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% */
1513 
1514     if (a <= sqrttwoe_6) {
1515                 Fdnest(z,5,&S6,&dS6dz);
1516         Fdnest(z,4,&S5,&dS5dz);
1517         Fdnest(z,3,&S4,&dS4dz);
1518         Fdnest(z,2,&S3,&dS3dz);
1519         Fdnest(z,1,&S2,&dS2dz);
1520         Fdnest(z,0,&S1,&dS1dz);
1521         dS1da = dS2dz-dS1dz;
1522         dS2da = dS3dz-dS2dz;
1523         dS3da = dS4dz-dS3dz;
1524         dS4da = dS5dz-dS4dz;
1525                 dS5da = dS6dz-dS5dz;
1526     }
1527 
1528     /* %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1529     %% Series expansion in a
1530     %% P.A. Robinson, J. Math. Phys. vol. 27, N�5, May 1986
1531     %% formula (10)
1532     %% valid for small a
1533     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% */
1534 
1535     if ((a>sqrttwoe_6)&(a <= 9.0)&(fabs(z-a)>6.0)&(fabs(z-a)<10*limit)) {
1536         if (a<=1.0) { /* some optimization for rel. error < 1e-12 */
1537             nmax = 13;
1538         } else if (a<=3) {
1539             nmax = 22;
1540         } else {
1541             nmax = 37;
1542         }
1543         z_a = z-a;
1544 
1545                 Fdnest(z_a,nmax+5,&F6,&dF6);
1546         Fdnest(z_a,nmax+4,&F5,&dF5);
1547         Fdnest(z_a,nmax+3,&F4,&dF4);
1548         Fdnest(z_a,nmax+2,&F3,&dF3);
1549         Fdnest(z_a,nmax+1,&F2,&dF2);
1550         Fdnest(z_a,nmax  ,&F1,&dF1);
1551 
1552         /* S1 = 0.0; */
1553         S2 = 0.0;
1554         S3 = 0.0;
1555         S4 = 0.0; /* S5 = 0.0; S6 = 0.0; */
1556         /* dS1dz = 0.0; */
1557         dS2dz = 0.0;
1558         dS3dz = 0.0;
1559         dS4dz = 0.0; /* dS5dz = 0.0; dS6dz = 0.0; */ 
1560         dS1da = 0.0;
1561         dS2da = 0.0;
1562         dS3da = 0.0;
1563         dS4da = 0.0;
1564         dS5da = 0.0;
1565 
1566         for (n = nmax;n>0;n--) {
1567             coeff = a/n;
1568             /* S1 = coeff*(S1+F1); */
1569             S2 = coeff*(S2+F2);
1570             S3 = coeff*(S3+F3);
1571                         S4 = coeff*(S4+F4);
1572                         /* S5 = coeff*(S5+F5);
1573                      S6 = coeff*(S6+F6); */
1574             /* dS1dz = coeff*(dS1dz+dF1) */;
1575             dS2dz = coeff*(dS2dz+dF2);
1576             dS3dz = coeff*(dS3dz+dF3);
1577             dS4dz = coeff*(dS4dz+dF4);
1578             /* dS5dz = coeff*(dS5dz+dF5);
1579                            dS6dz = coeff*(dS6dz+dF6); */
1580             dS1da = coeff*((dF2-dF1)+dS1da);
1581             dS2da = coeff*((dF3-dF2)+dS2da);
1582                 dS3da = coeff*((dF4-dF3)+dS3da);
1583             dS4da = coeff*((dF5-dF4)+dS4da);
1584                         dS5da = coeff*((dF6-dF5)+dS5da);
1585             /* F6 = F5;*/
1586                         F5 = F4;
1587             F4 = F3;
1588             F3 = F2;
1589             F2 = F1;
1590             dF6 = dF5;
1591             dF5 = dF4;
1592             dF4 = dF3;
1593             dF3 = dF2;
1594             dF2 = dF1;
1595             Fdnest(z_a,n-1,&F1,&dF1);
1596         }   /* end for*/
1597 
1598         expo = exp(-a);
1599         /* S1 = expo*(S1+F1); */
1600         S2 = expo*(S2+F2);
1601         S3 = expo*(S3+F3);
1602                 S4 = expo*(S4+F4);
1603                 /* S5 = expo*(S5+F5);
1604                S6 = expo*(S6+F6);*/
1605         /* dS1dz = expo*(dS1dz+dF1); */
1606         dS2dz = expo*(dS2dz+dF2);
1607         dS3dz = expo*(dS3dz+dF3);
1608         dS4dz = expo*(dS4dz+dF4);
1609         /* dS5dz = expo*(dS5dz+dF5); */
1610                 /* dS6dz = expo*(dS6dz+dF6); */
1611         dS1da = expo*(dS1da+(dF2-dF1));
1612         dS2da = expo*(dS2da+(dF3-dF2));
1613         dS3da = expo*(dS3da+(dF4-dF3));
1614         dS4da = expo*(dS4da+(dF5-dF4));
1615                 dS5da = expo*(dS5da+(dF6-dF5));
1616 
1617     } /* end if ((a>0.0)&(a<=9.0)&(fabs(z-a)>6.0)&(fabs(z-a)<10*limit)) */
1618 
1619     if ((a>0.0)&(a<=9.0)&(fabs(z-a)<=6.0)) {
1620         if (a<=1.0) {/* some optimization for rel. error < 1e-12 */
1621             nmax = 13;
1622         } else if (a<=3) {
1623             nmax = 22;
1624         } else {
1625             nmax = 37;
1626         }
1627         z_a = z-a;
1628         coeff = 1.0;
1629         Fdnest(z_a,0,&F1,&dF1);
1630         F2 = (1.0-z_a*F1)*0.66666666666666666667; /* use recursion relation */
1631         F3 = (1.0-z_a*F2)*0.4;
1632         F4 = (1.0-z_a*F3)*0.2857142857142857176;
1633         F5 = (1.0-z_a*F4)*0.22222222222222222222;
1634                 F6 = (1.0-z_a*F5)*0.18181818181818181818;
1635 
1636         dF2 = F2-F1; /* use derivation relation */
1637         dF3 = F3-F2;
1638         dF4 = F4-F3;
1639         dF5 = F5-F4;
1640                 dF6 = F6-F5;
1641         d2F1 = dF2-dF1; /* derivation relation again */
1642         d2F2 = dF3-dF2;
1643         d2F3 = dF4-dF3;
1644         d2F4 = dF5-dF4;
1645                 d2F5 = dF6-dF5;
1646         /* S1 = F1; */
1647         S2 = F2; 
1648         S3 = F3; 
1649         S4 = F4; 
1650         /* S5 = F5;
1651                S6 = F6; */ 
1652         /* dS1dz = dF1; */
1653         dS2dz = dF2;
1654         dS3dz = dF3; 
1655         dS4dz = dF4; 
1656         /* dS5dz = dF5;
1657                dS6dz = dF6; */
1658         dS1da = d2F1;
1659         dS2da = d2F2;
1660         dS3da = d2F3;
1661         dS4da = d2F4;
1662                 dS5da = d2F5;
1663         for (n = 1;n<=nmax;n++) {
1664                 /* F1 = F2; */ 
1665             F2 = F3;
1666             F3 = F4;
1667             F4 = F5;
1668             F5 = F6;
1669                 dF1 = dF2;
1670             dF2 = dF3;
1671             dF3 = dF4;
1672             dF4 = dF5;
1673             dF5 = dF6;
1674                 d2F1 = d2F2;
1675             d2F2 = d2F3;
1676             d2F3 = d2F4;
1677             d2F4 = d2F5;
1678         
1679                 F6 = (1.0-z_a*F6)/(n+5.5);
1680                 dF6=F6-F5;
1681                 d2F5 = dF6-dF5;
1682                 coeff*=(a/n);
1683                 /* S1 += (coeff*F1); */
1684                 S2 += (coeff*F2);
1685             S3 += (coeff*F3); 
1686             S4 += (coeff*F4);
1687             /* S5 += (coeff*F5);
1688                            S6 += (coeff*F6); */
1689                 /* dS1dz += (coeff*dF1); */
1690                 dS2dz += (coeff*dF2);
1691             dS3dz += (coeff*dF3); 
1692             dS4dz += (coeff*dF4); 
1693             /* dS5dz += (coeff*dF5);
1694                       dS6dz += (coeff*dF6); */
1695                 dS1da += (coeff*d2F1); 
1696             dS2da += (coeff*d2F2);
1697                 dS3da += (coeff*d2F3);
1698             dS4da += (coeff*d2F4);
1699                         dS5da += (coeff*d2F5);
1700         }
1701         expo = exp(-a);
1702 
1703         /* S1 *= expo; */
1704         S2 *= expo;
1705         S3 *= expo;
1706         S4 = expo;
1707         /*S5 *= expo;
1708                S6 *= expo;
1709         dS1dz *= expo; */
1710         dS2dz *= expo;
1711         dS3dz *= expo;
1712         dS4dz *= expo;
1713         /*   dS5dz *= expo;
1714                   dS6dz *= expo; */     
1715         dS1da *= expo;
1716         dS2da *= expo;
1717         dS3da *= expo;
1718         dS4da *= expo;
1719                 dS5da *= expo;
1720 
1721     } /* end if ((a>0.0)&(a<=9.0)&(fabs(z-a)<=6.0)) */  
1722 
1723 
1724     /* %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1725        %% Recursion formula
1726        %% P.A. Robinson, J. Math. Phys. vol. 27, N�5, May 1986
1727        %% formula (10)
1728        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% */
1729 
1730 
1731     if ((a>9.0)&(fabs(z-a)<=(10*limit))) {
1732         phi2 = a-z;
1733         psi = sqrt(a);
1734         /*   %%% Sign conventions for phi => 
1735              %%% See Brambilla p. 146 formula (16.51) */
1736    
1737         /*phi = conj(sqrt(phi2)); */
1738         if (phi2>0.0) {
1739             phi = sqrt(phi2);
1740             /* Z is real, contains no imaginary part */
1741             
1742 /*          Z1 = -2.0*dawson(psi-phi);*/
1743 /*          Z2 = -2.0*dawson(-psi-phi); */
1744 /*                        S0 = -0.5*(Z1+Z2)/phi;*/ /* F1_2 */
1745 /*          S1 = 0.5*(Z2-Z1)/psi;*/ /* F3_2 */
1746             dawson(psi-phi,&daw);
1747             rZ1 = -2.0*daw;
1748             dawson(-psi-phi,&daw);
1749             rZ2 = -2.0*daw; 
1750                         S0 = -0.5*(rZ1+rZ2)/phi; /* F1_2 */
1751             S1 = 0.5*(rZ2-rZ1)/psi; /* F3_2 */
1752 
1753                 } else if (phi2 == 0.0) {
1754                         S0 = 0.0;
1755             dawson(psi,&daw);
1756                         S1 = 2.0*daw/psi;
1757         } else { /* (phi2<0.0) */
1758             phi = -sqrt(-phi2);
1759 /*          Zfried (psi,-phi,&rZ1,&iZ1);*/
1760 /*          Zfried (-psi,-phi,&rZ2,&iZ2); */
1761 /*          only the real part is retained */
1762 /*          S0 = -0.5*(iZ1+iZ2)/phi;*/ /* F1_2 */
1763 /*          S1 = 0.5*(rZ2-rZ1)/psi;*/ /* F3_2 */
1764             ZFriedConte(psi-I*phi,&Z1);
1765             ZFriedConte(-psi-I*phi,&Z2);
1766             S0 = -0.5*(cimag(Z1)+cimag(Z2))/phi; /* F1_2 */
1767             S1 = 0.5*(creal(Z2)-creal(rZ1))/psi; /* F3_2 */
1768 
1769         }
1770   
1771         S2 = (1.0+phi2*S0-0.5*S1)/a; /* F_5/2 */
1772         S3 = (1.0+phi2*S1-1.5*S2)/a; /* F_7/2 */
1773         S4 = (1.0+phi2*S2-2.5*S3)/a; /* F_9/2 */
1774         S5 = (1.0+phi2*S3-3.5*S4)/a; /* F_11/2 */
1775                 S6 = (1.0+phi2*S4-4.5*S5)/a; /* F_13/2 */ 
1776         dS1dz = S1-S0;
1777         dS2dz = S2-S1;
1778         dS3dz = S3-S2;
1779         dS4dz = S4-S3;
1780         dS5dz = S5-S4;
1781                 dS6dz = S6-S5;
1782         dS1da = dS2dz-dS1dz;
1783         dS2da = dS3dz-dS2dz;
1784         dS3da = dS4dz-dS3dz;
1785         dS4da = dS5dz-dS4dz;
1786                 dS5da = dS6dz-dS5dz;
1787 
1788     } /* end if ((a>9.0)&(fabs(z-a)<=(10*limit))) */
1789 
1790 /* *RF3 = S1; */
1791     *RF5 = S2;
1792     *RF7 = S3;
1793     *RF9 = S4;
1794 /* *RF11 = S5; */
1795 /* *RF13 = S6; */
1796 /* *dRF3dz = dS1dz; */
1797     *dRF5dz = dS2dz;
1798     *dRF7dz = dS3dz;
1799     *dRF9dz = dS4dz;
1800 /* *dRF11dz = dS5dz; */
1801 /* *dRF13dz = dS6dz; */
1802     *dRF3da  = dS1da;
1803     *dRF5da  = dS2da;
1804     *dRF7da  = dS3da;
1805     *dRF9da  = dS4da;
1806     *dRF11da = dS5da;
1807 
1808 } 
1809 
1810 
1811 
1812 
1813 
1814 
1815 
1816 void ImShka(double z, double a, 
1817              double * IF5, double * IF7, double * IF9,
1818          double * dIF5dz, double * dIF7dz, double * dIF9dz,
1819              double * dIF3da, double * dIF5da, double * dIF7da, 
1820              double * dIF9da, double * dIF11da) {
1821 
1822 double IS1, IS3, IS5, IS7, IS9, IS11, IS13;
1823 double dIS3dz, dIS5dz, dIS7dz, dIS9dz, dIS11dz, dIS13dz;
1824 double dIS3da, dIS5da, dIS7da, dIS9da, dIS11da;
1825 double absz, sqrtone_za, argum, arg1, arg2, bes0, bes1, bes2;
1826 static double twosqrtpi = 3.5449077018110318818;
1827 static double sqrttwoe_6 = 1.414e-6;
1828 
1829 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1830 %% Imaginary part of Shkarovsky function F_q+3/2(z,a)
1831 %% along with its first partial derivatives.
1832 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1833 %% P.A. Robinson, J. Math. Phys. vol. 27, N�5, May 1986
1834 %% Formula (42) and (43)
1835 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% */
1836 
1837    if (z>=a) {
1838    /* No absorption */
1839     /* IS1 = 0.0 ; IS3 = 0.0 */
1840     IS5 = 0.0; IS7 = 0.0;  IS9 = 0.0; /* IS11 = 0.0; IS13 = 0.0; */
1841     /* dIS3dz = 0.0 */
1842     dIS5dz = 0.0;
1843     dIS7dz = 0.0; 
1844     dIS9dz = 0.0; /* dIS11dz = 0.0; dIS13dz = 0.0*/ 
1845     dIS3da = 0.0;
1846     dIS5da = 0.0;
1847     dIS7da = 0.0;
1848     dIS9da = 0.0;
1849     dIS11da = 0.0;
1850    } else if ((z/* a==0 Shkarovski function => Dnestrovski function %% */
1852 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1853 %% P.A. Robinson, J. Math. Phys. vol. 27, N�5, May 1986
1854 %% Formula (42)
1855 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% */
1856       absz = -z;
1857       IS3 = -twosqrtpi*sqrt(absz)*exp(-absz);
1858       IS1 = IS3 *0.5/absz;
1859       IS5 = IS3 *absz*0.66666666666666666667;
1860       IS7 = IS5 *absz*0.4;
1861       IS9 = IS7 *absz*0.2857142857142857176;
1862       IS11 = IS9 *absz*0.22222222222222222222;
1863       IS13 = IS11*absz*0.18181818181818181818;
1864       dIS3dz = IS3 - IS1; /* P.A. Robinson formula (27) */
1865       dIS5dz = IS5 - IS3;
1866       dIS7dz = IS7 - IS5;
1867       dIS9dz = IS9 - IS7;
1868       dIS11dz = IS11 - IS9;
1869       dIS13dz = IS13 - IS11;
1870       dIS3da = dIS5dz - dIS3dz; /* P.A. Robinson formula (5) and (27) */
1871       dIS5da = dIS7dz - dIS5dz;
1872       dIS7da = dIS9dz - dIS7dz;
1873       dIS9da = dIS11dz - dIS9dz;
1874       dIS11da = dIS13dz - dIS11dz;
1875    } else if ((z < a) & (a > sqrttwoe_6)) {
1876 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1877 %% P.A. Robinson, J. Math. Phys. vol. 27, N�5, May 1986
1878 %% Formula (43)
1879 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% */
1880       sqrtone_za = sqrt(1.0-z/a);
1881       argum = 2.0*a*sqrtone_za; /* Argument of modified Bessel function */
1882       arg1 = 1.0/argum;
1883       if (argum<25.0) { /* argum<25.0 => complete formula */
1884          /* Normalization factor sqrt(2*argum/pi) for modified Bessel
1885           functions I_n+1/2 is directly incorporated here */
1886          IS1 = -twosqrtpi*exp(z-2.0*a)*sqrt(a);
1887          IS3 = IS1 *sqrtone_za;
1888          IS5 = IS3 *sqrtone_za;
1889          IS7 = IS5 *sqrtone_za;
1890          IS9 = IS7 *sqrtone_za;
1891          IS11 = IS9 *sqrtone_za;
1892          IS13 = IS11*sqrtone_za;
1893          bes0 = cosh(argum)*arg1; /* Normalized modified Bessel func. I_-1/2 */
1894          IS1 *= bes0; 
1895          bes1 = sinh(argum)*arg1; /* Normalized modified Bessel func. I_+1/2 */
1896          IS3 *= bes1;
1897          arg2 = 0.5*argum*argum;
1898          if (argum>4.0e-2) { /* Normalized modified Bessel func. I_3/2 */
1899             bes2 = bes0-bes1*arg1; /* Large argument => use recurrence relation */
1900          } else { 
1901             /* Small argument => use series expansion */
1902             /* Precision of 1e-12 on I_3/2 is achieved */
1903             bes2 = argum*0.33333333333333333*
1904             (1.0+0.2*arg2*
1905             (1.0+arg2*0.0714285714285714294));
1906          }
1907          IS5 *= bes2;
1908          bes0 = bes1;
1909      bes1 = bes2;
1910          if (argum>4.0e-1) { /* Normalized modified Bessel function I_5/2 */
1911             bes2 = bes0-3.0*bes1*arg1;
1912          } else {
1913             bes2 = argum*argum*0.06666666666666666*
1914             (1.0+arg2*0.1428571428571428588*
1915             (1.0+arg2*0.05555555555555555*
1916             (1.0+arg2*0.030303030303030303*
1917             (1.0+arg2*0.0192307692307692322))));
1918          }
1919          IS7 *= bes2;
1920          bes0 = bes1; bes1=bes2;
1921          if (argum>1.2) { /* Normalized modified Bessel function I_7/2 */
1922             bes2 = bes0-5.0*bes1*arg1;
1923          } else {
1924             bes2 = argum*argum*argum*0.0095238095238095254*
1925             (1.0+arg2*0.111111111111111111*
1926             (1.0+arg2*0.045454545454545454*
1927             (1.0+arg2*0.0256410256410256413*
1928             (1.0+arg2*0.016666666666666666*
1929             (1.0+arg2*0.0117647058823529404*
1930             (1.0+arg2*0.0087719298245614034))))));
1931          }
1932          IS9 *= bes2;
1933          bes0 = bes1;
1934      bes1 = bes2;
1935      if (argum>1.9) {
1936             bes2 = bes0-7.0*bes1*arg1;
1937          } else { /* Normalized modified Bessel function I_9/2 */
1938             bes2 = argum*argum*argum*argum*0.001058201058201058201*
1939             (1.0+arg2*0.090909090909090909*
1940             (1.0+arg2*0.0384615384615384655*
1941             (1.0+arg2*0.0222222222222222222*
1942             (1.0+arg2*0.01470588235294117694*
1943             (1.0+arg2*0.0105263157894736846*
1944             (1.0+arg2*0.0079365079365079364*
1945             (1.0+arg2*0.0062111801242236022)))))));
1946          }
1947          IS11 *= bes2;
1948          bes0 = bes1;
1949      bes1 = bes2;
1950      if (argum>3.4) {
1951             bes2 = bes0-9.0*bes1*arg1;
1952          } else { /* Normalized modified Bessel function I_11/2 */
1953             bes2 = argum*argum*argum*argum*argum*9.62000962000962e-5*
1954             (1.0+arg2*0.0769230769230769250*
1955             (1.0+arg2*0.033333333333333333333333*
1956             (1.0+arg2*0.01960784313725490206*
1957             (1.0+arg2*0.0131578947368421042*
1958             (1.0+arg2*0.0095238095238095254*
1959             (1.0+arg2*0.0072463768115942032*
1960             (1.0+arg2*0.0057142857142857142*
1961             (1.0+arg2*0.0046296296296296295*
1962             (1.0+arg2*0.0038314176245210725)))))))));
1963          }
1964          IS13 *= bes2;
1965 
1966       } else { 
1967 /* (argum>=25.0), use asymptotic expression for Shkarovski function */
1968          arg2 = sqrt(a-z)-sqrt(a);
1969          IS1 = -0.25*twosqrtpi*exp(-arg2*arg2)/sqrt(a-z);
1970          IS3 = IS1 *sqrtone_za;
1971          IS5 = IS3 *sqrtone_za;
1972          IS7 = IS5 *sqrtone_za;
1973          IS9 = IS7 *sqrtone_za;
1974          IS11 = IS9 *sqrtone_za;
1975          IS13 = IS11*sqrtone_za;
1976          IS5 *= (1.0-arg1);
1977          IS7 *= (1.0+3.0*arg1*(-1.0+arg1));
1978          IS9 *= (1.0+arg1*(-6.0+15.0*arg1*(1.0-arg1)));
1979          IS11 *= (1.0+arg1*(-10.0+arg1*(45.0+105.0*arg1*(arg1-1.0))));
1980          IS13 *= (1.0+arg1*(-15.0+arg1*(105.0+arg1*(-420.0+945.0*arg1*(1.0-arg1))))); 
1981       } /* end if (argum>=25) */
1982       dIS3dz = IS3 -IS1; /* P.A. Robinson formula (27) */
1983       dIS5dz = IS5 -IS3;
1984       dIS7dz = IS7 -IS5;
1985       dIS9dz = IS9 -IS7;
1986       dIS11dz = IS11-IS9;
1987       dIS13dz = IS13-IS11;
1988       dIS3da = dIS5dz -dIS3dz; /* P.A. Robinson formula (5) and (27) */
1989       dIS5da = dIS7dz -dIS5dz;
1990       dIS7da = dIS9dz -dIS7dz;
1991       dIS9da = dIS11dz-dIS9dz;
1992       dIS11da = dIS13dz-dIS11dz;
1993    } /* end if (a>sqrttwoe_6) &(z<a) */
1994 
1995    /* * IF1 = IS1;
1996    * IF3 = IS3; */   
1997    * IF5 = IS5;
1998    * IF7 = IS7;
1999    * IF9 = IS9;
2000    /* 
2001    * IF11 = IS11; 
2002    * IF13 = IS13;
2003    * dIF3dz = dIS3dz; */   
2004    * dIF5dz = dIS5dz;
2005    * dIF7dz = dIS7dz;
2006    * dIF9dz = dIS9dz;
2007    /* * dIF11dz = dIS11dz; 
2008     * dIF13dz = dIS13dz; */
2009    * dIF3da = dIS3da;
2010    * dIF5da = dIS5da;
2011    * dIF7da = dIS7da;
2012    * dIF9da = dIS9da;
2013    * dIF11da = dIS11da;
2014 
2015 }
2016 
2017 
2018 void Fdnest(z,q,F,dFdz)
2019 double z;
2020 int q;
2021 double *F,*dFdz;
2022 
2023 {
2024 /* Computes real part of Dnestrovskii's function F_q+3/2(z), for q>=0;
2025  along with its first derivative for q>=0 */
2026 
2027     double Fr,dFrdz,oldFr,ksi,multi,z_1,expo,lmt1,lmt3,daw;
2028     int n, nmx,nnmx,nnnmx;
2029     static double twosqrtpi = 3.54490770181103188818;
2030 
2031     static int taille = 46;
2032     static double limit1[46] = {-31.5,-34.6,-38.0,-39.7,-42.4,
2033                                 -44.8,-47.3,-49.4,-51.4,-53.4,
2034                                 -55.5,-57.5,-59.5,-61.5,-63.5,-65.5,-67.4,-68.7,-70.5,-72.6,-74.4,
2035                                 -76.3,-77.5,-79.7,-81.8,-82.7,-84.8,-86.4,-87.7,-89.5,-91.6,-92.9,
2036                                 -94.5,-96.6,-97.5,-99.6,-100.5,-102.6,-103.5,-105.5,-107.25,
2037                                 -108.6,-110.7,-111.6,-113.7,-114.6};
2038     static double limit3[46] = {40.0,14.0,9.0,7.0,6.0,6.0,6.0,6.0,
2039                                 6.2,6.0,7.0,7.0,6.5,7.5,7.5,
2040                                 7.0,8.5,8.5,8.0,8.0,9.5,9.0,9.0,8.5,11.2,11.2,11.0,11.0,10.0,10.0,
2041                                 13.5,13.5,13.2,13.0,13.0,12.5,12.0,12.0,11.5,11.0,10.7,17.0,17.0,16.5,
2042                                 16.0,16.0};
2043     static double lmt2=0.0;
2044     static int nmax[46] = {31,33,36,36,38,39,41,42,43,44,45,46,47,48,
2045                            49,50,51,51,52,53,54,55,55,56,57,57,58,59,59,60,61,61,62,63,63,64,64,65,
2046                            65,66,67,67,68,68,69,69};
2047     static int nnmax[46] = {78,83,87,90,95,98,102,105,108,111,114,117,120,
2048                             123,126,129,132,134, 136,139,142,144,146,149,152,153,156,158,160,162,165,
2049                             167,169,172,173,176,177,180,181,184,186,188,191,192,195,196};
2050     static int nnnmax[46] = {4,8,11,14,16,16,16,16,15,15,14,14,14,13,13,13,
2051                              12,12,12,12,11,11,11,11,10,10,10,10,10,10,9,9,9,9,9,9,9,9,9,9,9,8,8,8,8,8};
2052 
2053     if (q < taille) { /* optimized for relative error < 1e-12 on F */
2054         lmt1 = limit1[q+1];
2055         lmt3 = limit3[q+1];
2056         nmx = nmax[q+1];
2057         nnmx = nnmax[q+1];
2058         nnnmx = nnnmax[q+1];
2059     } else { 
2060 /* not optimized, relative error of 1e-12 garantied for z<0 */
2061         lmt1 = -43.3891 - 1.5934*q - 5.0; /* apply a margin of 5 */
2062         lmt3 = limit3[taille];            
2063         nmx = (int)(22.1786+7.0234*sqrt(q));
2064         nnmx = (int) (100.0584 + 2.1494*q) + 10; /* apply a margin of 10 */
2065         nnnmx = nnnmax[taille];
2066     }
2067 
2068     if (z/*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2070 %% Asymptotic approximation for |z+q|>>1
2071 %% P.A. Robinson, J. Math. Phys. vol. 27, N�5, May 1986
2072 %% Formula (32)
2073 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2074 */
2075         z_1=1.0/z;
2076         Fr=0.0;
2077         dFrdz=0.0;
2078         multi=nmx+q-0.5;
2079         for (n = nmx;n > 1;n--) {
2080         Fr = -z_1*multi*(1.0 + Fr);
2081             dFrdz = z_1*multi*(n*z_1 - dFrdz);
2082             multi--;
2083         }
2084     Fr = z_1*(1.0 + Fr);
2085     dFrdz = z_1*(dFrdz-z_1);
2086 
2087     } else if ((z >= lmt1) & (z < lmt2)) {
2088 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2089 %% Series development (unstable for positive z ?)
2090 %% P.A. Robinson, J. Math. Phys. vol. 27, N�5, May 1986
2091 %% Formula (30) for z<0
2092 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2093 */
2094         z_1 = -z;
2095         Fr=0.0;
2096         dFrdz=0.0;
2097         multi=nnmx-0.5-q;
2098         for (n=nnmx;n>0;n--) {
2099             Fr=z_1/n*(1.0/multi+Fr);
2100             dFrdz=z_1/n*(1.0/multi/(multi+1.0)+dFrdz);
2101             multi--;
2102         }
2103         expo=exp(z);
2104 /* !!!! Formula only valid for z<=0 
2105    !!!! Extra terms should be added for z>0 !!!! */
2106     
2107         Fr=-expo*(Fr+1.0/multi);
2108         dFrdz=-expo*(dFrdz+1.0/multi/(multi+1.0));
2109 
2110     } else if ((z >= lmt2) & (z < lmt3)) {
2111 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2112 %% Recurrence relation (unstable for large q and z)
2113 %% P.A. Robinson, J. Math. Phys. vol. 27, N�5, May 1986
2114 %% Formula (19)
2115 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2116 */   
2117 /* First compute F_3/2 */
2118 
2119         if (z == 0.0) {
2120             ksi = sqrt(z);
2121             Fr = 2.0-twosqrtpi*ksi*erfc(ksi)*exp(z);
2122     } else {
2123             ksi = -sqrt(-z);
2124         dawson(ksi,&daw);
2125             Fr = 2.0 - 4.0*ksi*daw;
2126         }
2127         
2128     /* Then apply recurrence relation F_q+1/2 => F_q+3/2
2129    ! recurrence relation becomes unstable for large |z| arguments ! */
2130         oldFr = (1.0 - 0.5*Fr)/(z + 1.0e-16);
2131         for (n = 1;n <= q;n++) {
2132         if (n == q) {
2133             oldFr=Fr;
2134             }     
2135             Fr = (1.0 - z*Fr)/(n + 0.5);
2136         }
2137         dFrdz=Fr-oldFr;
2138 
2139     } else { /*if (z>=lmt3) */
2140 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2141 %% Continued fraction (unstable for negative z)
2142 %% P.A. Robinson, J. Math. Phys. vol. 27, N�5, May 1986
2143 %% Formula (34)
2144 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2145 */
2146         multi = q + 1.5 + nnnmx;
2147         Fr = 0.0;
2148         oldFr = 0.0;
2149         for (n = nnnmx;n > 0;n--) {
2150             Fr = n/(z + multi/(1.0 + Fr));
2151             multi--;
2152             oldFr = n/(z + multi/(1.0 + oldFr));
2153         }
2154         Fr = 1.0/(z + multi/(1.0 + Fr));
2155         multi--;
2156         dFrdz = Fr - 1.0/(z + multi/(1.0 + oldFr));  
2157 
2158     } /* end if z */
2159 
2160     *F = Fr;
2161     *dFdz = dFrdz;   
2162    
2163 }
2164 
2165 
2166 
2167 void ImFdnest(z,q,ImF,dImF)
2168 double z;
2169 int q; 
2170 double *ImF,*dImF;
2171 {
2172 double IFdnest,dIFdz,absz,oldIF;
2173 int n;
2174 static double twosqrtpi = 3.5449077018110318818;
2175 
2176 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2177 % [IF,dIFdz]=ImFdnest(z,q);
2178 % Imaginary part of Dnestrovskii's function F_q+3/2(z), 
2179 % for integer q>=0;
2180 % along with its first derivative for q>=1.
2181 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2182 %% P.A. Robinson, J. Math. Phys. vol. 27, N�5, May 1986
2183 %% Formula (42) for z<0 plus recursion
2184 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2185 */
2186 
2187     if (z < 0.0) {
2188         absz=-z;
2189         IFdnest = -twosqrtpi*sqrt(absz)*exp(-absz);
2190         oldIF = IFdnest*0.5/absz;
2191         for (n = 1;n <= q;n++) {
2192             if (n == q) {
2193                 oldIF = IFdnest;
2194             }
2195             IFdnest *= (absz/(n+0.5));
2196         }
2197         dIFdz = IFdnest - oldIF;
2198     } else {
2199         IFdnest = 0.0;
2200         dIFdz = 0.0;
2201    }
2202    *ImF = IFdnest;
2203    *dImF = dIFdz;
2204 } 
2205 
2206 
2207