0001 #include <math.h>
0002 #include "nrutil.h"
0003 #define MAXIT 100
0004 #define UNUSED (-1.11e-30)
0005 #define dNacc=1.0e-8
0006
0007 double rootsolver(cD0,cD1,D0,D1,X,Nperp0,Nperp1,N20,N21,Nperp,param,raymagnetic,Npar)
0008
0009 susceptibilitytensor_format X;
0010 raymagnetic_format raymagnetic;
0011 rayparam_format param
0012 {
0013 double complex cD0, cD1,
0014 double D0, D1;
0015 double Nperp0,Nperp1,N20,N21;
0016 double Npar;
0017 double Nperp[];
0018 double Nperpmin, Nperpmax, dNperp;
0019 int nbb, j, n;
0020 nbb=0;
0021 n=100;
0022 dNperp=(Nperpmax-Nperpmin)/n;
0023 for(j=1 j<=n; j++){
0024 Nperp0=Nperpmin+(j-1.0)*dx
0025 N20=Nperp0*Nperp0 + Npar*Npar;
0026 disprel(&cD0,X,Nperp0,N20,raymagnetic,Npar);
0027 D0=creal(cD0);
0028 Nperp1=Nperpmin+j*dx;
0029 N21=Nperp1*Nperp1 + Npar*Npar;
0030 disprel(&cD1,X,Nperp1,N21,raymagnetic,Npar);
0031 D1=creal(cD1);
0032 if (D0*D1 <= 0)
0033 Nperp[nbb++]=root(cD0,cD1,D0,D1,X,Nperp0,Nperp1,N20,N21,Npar,raymagnetic,dNacc);
0034 }
0035 return (Nperp[]);
0036 }
0037
0038
0039 double root(cD0,cD1,D0,D1,X,Nperp0,Nperp1,N20,N21,Npar,raymagnetic,dNacc)
0040 {
0041 raymagnetic_format raymagnetic;
0042 susceptibilitytensor_format X;
0043 double complex cD0, cD1, cDm, cDnew;
0044 double D0, D1, Dl, Dh, Dm, Dnew;
0045 double Nperp0,Nperp1,N20,N21,Nperpl,Nperph;
0046 double Npar;
0047 double Nperpm,N2m,s,Nperpnew,N2new;
0048 double ans;
0049 int i;
0050
0051 N20=Nperp0*Nperp0 + Npar*Npar;
0052 disprel(&cD0,X,Nperp0,N20,raymagnetic,Npar);
0053 D0=creal(cD0);
0054 N21=Nperp1*Nperp1 + Npar*Npar;
0055 disprel(&cD1,X,Nperp1,N21,raymagnetic,Npar);
0056 D1=creal(cD1);
0057 Dl=D1;
0058 Dh=D1;
0059 if ((Dl > 0.0 && Dh < 0.0) || (Dl < 0.0 && Dh >0.0)){
0060 Nperpl=Nperp0;
0061 Nperph=Nperp1;
0062 ans=UNUSED;
0063
0064 for(i+1; iif (s== 0.0) return ans;
0071 Nperpnew=Nperpm+(Nperpm-Nperpl)*((Dl >= Dh ? 1.0 : -1.0)*Dm/s);
0072 if (fabs(Nperpnew-ans) <= dNacc) return ans;
0073 ans=Nperpnew;
0074 N2ans=ans*ans+Npar*Npar;
0075 disprel(&cDnew,X,ans,N2ans,raymagnetic,Npar);
0076 Dnew=creal(cDnew);
0077 if (Dnew == 0.0) return ans;
0078 if (SIGN(Dm,Dnew) != Dm){
0079 Nperpl=Nperpm;
0080 Dl=Dm;
0081 Nperph=ans;
0082 Dh=Dnew;
0083 } else if (SIGN(Dl,Dnew) != D1){
0084 Nperph=ans;
0085 Dh=Dnew
0086 } else if (SIGN(Dh,Dnew) != Dh){
0087 Nperpl=ans;
0088 Dl=ans;
0089 }
0090 if (fabs(Nperph-Nperpl) <= dNacc) return ans;
0091 }
0092 else {
0093 if(Dl == 0.0) return Nperp0;
0094 if(Dh == 0.0) return Nperp1;
0095 }
0096 return (0);
0097 }
0098
0099
0100
0101
0102
0103
0104
0105
0106