55static const int NSD = 7;
143 for(
int j=0; j <
NAX; j++ )
150 for(
int j=0; j <
NDAT; j++ )
155 vector< complex<double> >
n[
NAX];
211 double*,
double*,
double*,
int*),
212 double*,
double*,
double*,
int*);
217 double*,
double*,
double*,
int*);
230inline double Drude(
double,
double,
double,
double);
239 long,
long,
int,
bool,
bool*);
242STATIC void init_eps(
double,
long,
const vector<grain_data>&,vector< complex<double> >&);
244 void(*)(complex<double>,
const vector<double>&,
const vector< complex<double> >&,
245 long,complex<double>*,
double*,
double*),
246 const vector<double>&,
const vector< complex<double> >&,
long,complex<double>,
double);
247STATIC void Stognienko(complex<double>,
const vector<double>&,
const vector< complex<double> >&,
248 long,complex<double>*,
double*,
double*);
249STATIC void Bruggeman(complex<double>,
const vector<double>&,
const vector< complex<double> >&,
250 long,complex<double>*,
double*,
double*);
265 double*,
double*,
long*);
266STATIC void anomal(
double,
double*,
double*,
double*,
double*,
double,
double);
267STATIC void bigk(complex<double>,complex<double>*);
273 const char *szd_file,
302 const long NPTS_TABLE = 100L;
313 fprintf(
ioQQQ,
" Illegal number of size distribution bins: %ld\n", sd.
nPart );
314 fprintf(
ioQQQ,
" The number should be between 1 and 99.\n" );
319 vector<grain_data> gdArr(2);
323 string rfi_string ( rfi_file );
324 if( rfi_string.find(
".rfi" ) != string::npos )
328 else if( rfi_string.find(
".mix" ) != string::npos )
334 fprintf(
ioQQQ,
" Refractive index file name %s has wrong extention\n", rfi_file );
335 fprintf(
ioQQQ,
" It should have extention .rfi or .mix.\n" );
341 fprintf(
ioQQQ,
" Illegal number of size distribution bins: %ld\n", sd.
nPart );
342 fprintf(
ioQQQ,
" The number should always be 1 for OPC_TABLE files.\n" );
347 fprintf(
ioQQQ,
" Illegal value for the specific density: %.4e\n", gd.
rho );
352 fprintf(
ioQQQ,
" Illegal value for the molecular weight: %.4e\n", gd.
mol_weight );
359 strcpy(chFile,rfi_file);
365 strcpy(chFile2,szd_file);
373 sprintf(ext,
"%02ld",nbin);
374 strcat(strcat(strcat(strcat(strcat(chFile,
"_"),chFile2),
"_"),ext),
".opc");
378 strcat(strcat(strcat(chFile,
"_"),chFile2),
".opc");
391 vector<double> inv_att_len(
rfield.nupper );
393 fprintf(
ioQQQ,
"\n Starting mie_write_opc, output will be written to %s\n\n", chFile );
399 lgErr = lgErr || ( fprintf(fdes,
"# this file was created by Cloudy %s (%s) on %s",
401 lgErr = lgErr || ( fprintf(fdes,
"# ===========================================\n#\n") < 0 );
402 lgErr = lgErr || ( fprintf(fdes,
"%12ld # magic number opacity file\n",
MAGIC_OPC) < 0 );
403 lgErr = lgErr || ( fprintf(fdes,
"%12ld # magic number rfi/mix file\n",gd.
magic) < 0 );
404 lgErr = lgErr || ( fprintf(fdes,
"%12ld # magic number szd file\n",sd.
magic) < 0 );
408 strncpy(hlp1,chFile,(
size_t)(
LABELSUB1+1));
415 strncpy(hlp2,chFile2,(
size_t)(
LABELSUB2+1));
422 strcpy(chGrainLabel,
" ");
427 strcat(strcat(strcat(strcat(chGrainLabel,hlp1),
"-"),hlp2),
"xx");
428 lgErr = lgErr || ( fprintf(fdes,
"%-12.12s # grain type label, xx will be replaced by bin no.\n",
433 strcat(strcat(strcat(chGrainLabel,hlp1),
"-"),hlp2);
434 lgErr = lgErr || ( fprintf(fdes,
"%-12.12s # grain type label\n", chGrainLabel) < 0 );
437 lgErr = lgErr || ( fprintf(fdes,
"%.6e # specific weight (g/cm^3)\n",gd.
rho) < 0 );
438 lgErr = lgErr || ( fprintf(fdes,
"%.6e # molecular weight of grain molecule (amu)\n",gd.
mol_weight) < 0 );
439 lgErr = lgErr || ( fprintf(fdes,
"%.6e # average molecular weight per atom (amu)\n", gd.
atom_weight) < 0 );
440 lgErr = lgErr || ( fprintf(fdes,
"%.6e # abundance of grain molecule at max depletion\n",gd.
abun) < 0 );
441 lgErr = lgErr || ( fprintf(fdes,
"%.6e # depletion efficiency\n",gd.
depl) < 0 );
442 lgErr = lgErr || ( fprintf(fdes,
"%.6e # average grain radius <a^3>/<a^2>, full size distr (cm)\n",
444 lgErr = lgErr || ( fprintf(fdes,
"%.6e # average grain surface area <4pi*a^2>, full size distr (cm^2)\n",
446 lgErr = lgErr || ( fprintf(fdes,
"%.6e # average grain volume <4/3pi*a^3>, full size distr (cm^3)\n",
448 lgErr = lgErr || ( fprintf(fdes,
"%.6e # total grain radius Int(a) per H, full size distr (cm/H)\n",
450 lgErr = lgErr || ( fprintf(fdes,
"%.6e # total grain area Int(4pi*a^2) per H, full size distr (cm^2/H)\n",
452 lgErr = lgErr || ( fprintf(fdes,
"%.6e # total grain vol Int(4/3pi*a^3) per H, full size distr (cm^3/H)\n",
454 lgErr = lgErr || ( fprintf(fdes,
"%.6e # work function (Ryd)\n",gd.
work) < 0 );
455 lgErr = lgErr || ( fprintf(fdes,
"%.6e # gap between valence and conduction band (Ryd)\n",gd.
bandgap) < 0 );
456 lgErr = lgErr || ( fprintf(fdes,
"%.6e # efficiency of thermionic emission\n",gd.
therm_eff) < 0 );
457 lgErr = lgErr || ( fprintf(fdes,
"%.6e # sublimation temperature (K)\n",gd.
subl_temp) < 0 );
458 lgErr = lgErr || ( fprintf(fdes,
"%12d # material type, 1=carbonaceous, 2=silicate\n",gd.
matType) < 0 );
459 lgErr = lgErr || ( fprintf(fdes,
"#\n# abundances of constituent elements rel. to hydrogen\n#\n") < 0 );
461 for( nelem=0; nelem <
LIMELM; nelem++ )
463 lgErr = lgErr || ( fprintf(fdes,
"%.6e # %s\n",gd.
elmAbun[nelem],
469 lgErr = lgErr || ( fprintf(fdes,
"#\n# entire size distribution, amin=%.5f amax=%.5f micron\n",
471 lgErr = lgErr || ( fprintf(fdes,
"#\n%.6e # ratio a_max/a_min in each size bin\n",
473 lgErr = lgErr || ( fprintf(fdes,
"#\n# size distribution function\n#\n") < 0 );
474 lgErr = lgErr || ( fprintf(fdes,
"%12ld # number of table entries\n#\n",NPTS_TABLE+1) < 0 );
475 lgErr = lgErr || ( fprintf(fdes,
"# ============================\n") < 0 );
476 lgErr = lgErr || ( fprintf(fdes,
"# size (micr) a^4*dN/da (cm^3/H)\n#\n") < 0 );
477 for( i=0; i <= NPTS_TABLE; i++ )
483 lgErr = lgErr || ( fprintf(fdes,
"%.6e %.6e\n",
radius,a4dNda) < 0 );
488 lgErr = lgErr || ( fprintf(fdes,
"#\n") < 0 );
489 lgErr = lgErr || ( fprintf(fdes,
"%.6e # a_max/a_min = 1 for single sized grain\n", 1. ) < 0 );
490 lgErr = lgErr || ( fprintf(fdes,
"%12ld # no size distribution table\n",0L) < 0 );
498 lgErr = lgErr || ( fprintf(fdes,
"#\n") < 0 );
499 lgErr = lgErr || ( fprintf(fdes,
"%s # check 1\n",
continuum.mesh_md5sum.c_str()) < 0 );
501 if(
cpu.i().big_endian() )
502 lgErr = lgErr || ( fprintf(fdes,
"%23.8x %8.8x # check 2\n",u.i[0],u.i[1]) < 0 );
504 lgErr = lgErr || ( fprintf(fdes,
"%23.8x %8.8x # check 2\n",u.i[1],u.i[0]) < 0 );
505 u.x = double(
rfield.egamry);
506 if(
cpu.i().big_endian() )
507 lgErr = lgErr || ( fprintf(fdes,
"%23.8x %8.8x # check 3\n",u.i[0],u.i[1]) < 0 );
509 lgErr = lgErr || ( fprintf(fdes,
"%23.8x %8.8x # check 3\n",u.i[1],u.i[0]) < 0 );
511 if(
cpu.i().big_endian() )
512 lgErr = lgErr || ( fprintf(fdes,
"%23.8x %8.8x # check 4\n",u.i[0],u.i[1]) < 0 );
514 lgErr = lgErr || ( fprintf(fdes,
"%23.8x %8.8x # check 4\n",u.i[1],u.i[0]) < 0 );
515 lgErr = lgErr || ( fprintf(fdes,
"%32ld # rfield.nupper\n",
rfield.nupper) < 0 );
516 lgErr = lgErr || ( fprintf(fdes,
"%32ld # number of size distr. bins\n#\n",sd.
nPart) < 0 );
530 vector<int> ErrorIndex(
rfield.nupper );
532 for( p=0; p < sd.
nPart; p++ )
546 volfrac = sd.
vol*frac/volnorm;
547 fprintf(
ioQQQ,
" Starting size bin %ld, amin=%.5f amax=%.5f micron\n",
549 lgErr = lgErr || ( fprintf(fdes,
"# size bin %ld, amin=%.5f amax=%.5f micron\n",
551 lgErr = lgErr || ( fprintf(fdes,
"%.6e # average grain ",3.*sd.
vol/sd.
area) < 0 );
552 lgErr = lgErr || ( fprintf(fdes,
"radius <a^3>/<a^2>, this bin (cm)\n") < 0 );
553 lgErr = lgErr || ( fprintf(fdes,
"%.6e # average ",sd.
area) < 0 );
554 lgErr = lgErr || ( fprintf(fdes,
"grain area <4pi*a^2>, this bin (cm^2)\n") < 0 );
555 lgErr = lgErr || ( fprintf(fdes,
"%.6e # average ",sd.
vol) < 0 );
556 lgErr = lgErr || ( fprintf(fdes,
"grain volume <4/3pi*a^3>, this bin (cm^3)\n") < 0 );
557 lgErr = lgErr || ( fprintf(fdes,
"%.6e # total grain ",sd.
radius*frac/gd.
norm) < 0 );
558 lgErr = lgErr || ( fprintf(fdes,
"radius Int(a) per H, this bin (cm/H)\n") < 0 );
559 lgErr = lgErr || ( fprintf(fdes,
"%.6e # total grain area ",sd.
area*frac/gd.
norm) < 0 );
560 lgErr = lgErr || ( fprintf(fdes,
"Int(4pi*a^2) per H, this bin (cm^2/H)\n") < 0 );
561 lgErr = lgErr || ( fprintf(fdes,
"%.6e # total grain volume ",sd.
vol*frac/gd.
norm) < 0 );
562 lgErr = lgErr || ( fprintf(fdes,
"Int(4/3pi*a^3) per H, this bin (cm^3/H)\n#\n") < 0 );
565 lgErrorOccurred =
false;
568 for( i=0; i <
rfield.nupper; i++ )
583 ErrorIndex[i] =
max(ErrorIndex[i],Error);
584 acs_abs[p][i] += cs_abs*gd.
wt[gd.
cAxis];
585 acs_sct[p][i] += cs_sct*gd.
wt[gd.
cAxis];
586 a1g[p][i] += cs_sct*(1.-cosb)*gd.
wt[gd.
cAxis];
588 lgErrorOccurred =
mie_auxiliary2(ErrorIndex,acs_abs,acs_sct,a1g,p,i);
593 ErrorIndex[i] =
min(Error,2);
594 lgErrorOccurred = lgErrorOccurred || ( Error > 0 );
595 acs_abs[p][i] = cs_abs*gd.
norm;
596 acs_sct[p][i] = cs_sct*gd.
norm;
601 acs_abs[p][i] = 1.3121e-23*volfrac*gd.
norm;
602 acs_sct[p][i] = 2.6242e-23*volfrac*gd.
norm;
610 ErrorIndex[i] =
max(ErrorIndex[i],Error);
611 acs_abs[p][i] += cs_abs*gd2.
wt[gd2.
cAxis];
612 acs_sct[p][i] += 0.1*cs_abs*gd2.
wt[gd2.
cAxis];
613 a1g[p][i] += 0.1*cs_abs*1.*gd2.
wt[gd2.
cAxis];
615 lgErrorOccurred =
mie_auxiliary2(ErrorIndex,acs_abs,acs_sct,a1g,p,i);
625 ErrorIndex[i] =
max(ErrorIndex[i],Error);
626 acs_abs[p][i] += cs_abs*gd2.
wt[gd2.
cAxis];
627 acs_sct[p][i] += 0.1*cs_abs*gd2.
wt[gd2.
cAxis];
628 a1g[p][i] += 0.1*cs_abs*1.*gd2.
wt[gd2.
cAxis];
630 lgErrorOccurred =
mie_auxiliary2(ErrorIndex,acs_abs,acs_sct,a1g,p,i);
640 ErrorIndex[i] =
max(ErrorIndex[i],Error);
641 acs_abs[p][i] += cs_abs*gd2.
wt[gd2.
cAxis];
642 acs_sct[p][i] += 0.1*cs_abs*gd2.
wt[gd2.
cAxis];
643 a1g[p][i] += 0.1*cs_abs*1.*gd2.
wt[gd2.
cAxis];
645 lgErrorOccurred =
mie_auxiliary2(ErrorIndex,acs_abs,acs_sct,a1g,p,i);
648 fprintf(
ioQQQ,
" Insanity in mie_write_opc\n" );
655 if( lgErrorOccurred )
657 strcpy(chString,
"absorption cs");
659 strcpy(chString,
"scattering cs");
661 strcpy(chString,
"asymmetry parameter");
665 for( i=0; i <
rfield.nupper; i++ )
667 acs_abs[p][i] /= gd.
norm;
670 acs_sct[p][i] /= gd.
norm;
674 lgErr = lgErr || ( fprintf(fdes,
"#\n") < 0 );
675 lgErr = lgErr || ( fprintf(fdes,
"# ===========================================\n") < 0 );
676 lgErr = lgErr || ( fprintf(fdes,
"# anu (Ryd) abs_cs_01 (cm^2/H) abs_cs_02.....\n#\n") < 0 );
678 for( i=0; i <
rfield.nupper; i++ )
680 lgErr = lgErr || ( fprintf(fdes,
"%.6e ",
rfield.anu[i]) < 0 );
681 for( p=0; p < sd.
nPart; p++ )
683 lgErr = lgErr || ( fprintf(fdes,
"%.6e ",acs_abs[p][i]) < 0 );
685 lgErr = lgErr || ( fprintf(fdes,
"\n") < 0 );
688 lgErr = lgErr || ( fprintf(fdes,
"#\n") < 0 );
689 lgErr = lgErr || ( fprintf(fdes,
"# ===========================================\n") < 0 );
690 lgErr = lgErr || ( fprintf(fdes,
"# anu (Ryd) sct_cs_01 (cm^2/H) sct_cs_02.....\n#\n") < 0 );
692 for( i=0; i <
rfield.nupper; i++ )
694 lgErr = lgErr || ( fprintf(fdes,
"%.6e ",
rfield.anu[i]) < 0 );
695 for( p=0; p < sd.
nPart; p++ )
697 lgErr = lgErr || ( fprintf(fdes,
"%.6e ",acs_sct[p][i]) < 0 );
699 lgErr = lgErr || ( fprintf(fdes,
"\n") < 0 );
702 lgErr = lgErr || ( fprintf(fdes,
"#\n") < 0 );
703 lgErr = lgErr || ( fprintf(fdes,
"# ===========================================\n") < 0 );
704 lgErr = lgErr || ( fprintf(fdes,
"# anu (Ryd) (1-g)_bin_01 (1-g)_bin_02.....\n#\n") < 0 );
706 for( i=0; i <
rfield.nupper; i++ )
708 lgErr = lgErr || ( fprintf(fdes,
"%.6e ",
rfield.anu[i]) < 0 );
709 for( p=0; p < sd.
nPart; p++ )
712 lgErr = lgErr || ( fprintf(fdes,
"%.6e ",
min(a1g[p][i],1.)) < 0 );
714 lgErr = lgErr || ( fprintf(fdes,
"\n") < 0 );
717 fprintf(
ioQQQ,
" Starting calculation of inverse attenuation length\n" );
718 strcpy(chString,
"inverse attenuation length");
735 fprintf(
ioQQQ,
" mie_write_opc detected unknown ial type: %d\n" , icase );
744 lgErr = lgErr || ( fprintf(fdes,
"#\n") < 0 );
745 lgErr = lgErr || ( fprintf(fdes,
"# ===========================================\n") < 0 );
746 lgErr = lgErr || ( fprintf(fdes,
"# anu (Ryd) inverse attenuation length (cm^-1)\n#\n") < 0 );
748 for( i=0; i <
rfield.nupper; i++ )
750 lgErr = lgErr || ( fprintf(fdes,
"%.6e %.6e\n",
rfield.anu[i],inv_att_len[i]) < 0 );
757 fprintf(
ioQQQ,
"\n Error writing file: %s\n", chFile );
758 if( remove(chFile) == 0 )
760 fprintf(
ioQQQ,
" The file has been removed\n" );
766 fprintf(
ioQQQ,
"\n Opacity file %s written succesfully\n\n", chFile );
769 fprintf(
ioQQQ,
"\n !!! Warnings were detected !!!\n\n" );
787 const double TOLER = 1.e-3;
790 if( strcmp(auxCase,
"init") == 0 )
792 double mass,
radius, CpMolecule;
809 fprintf(
ioQQQ,
"\n This size distribution can only be combined with"
810 " carbonaceous material, bailing out...\n" );
844 delta = fabs(sd->
vol-oldvol)/sd->
vol;
846 }
while( sd->
nmul <= 1024 && delta >
TOLER );
850 fprintf(
ioQQQ,
" could not converge integration of size distribution\n" );
861 fprintf(
ioQQQ,
" insane case for grain size distribution: %d\n" , sd->
sdCase );
866 else if( strcmp(auxCase,
"step") == 0 )
883 step = (amax - amin)/(
double)sd->
nPart;
884 amin = amin + (double)sd->
cPart*step;
885 amax =
min(amax,amin + step);
898 fprintf(
ioQQQ,
" insane case for grain size distribution: %d\n" , sd->
sdCase );
905 fprintf(
ioQQQ,
" mie_auxiliary called with insane argument: %s\n", auxCase );
922 bool lgErrorOccurred =
false;
923 if( ErrorIndex[i] > 0 )
925 ErrorIndex[i] =
min(ErrorIndex[i],2);
926 lgErrorOccurred =
true;
929 switch( ErrorIndex[i] )
942 a1g[p][i] /= acs_sct[p][i];
945 fprintf(
ioQQQ,
" Insane value for ErrorIndex: %d\n", ErrorIndex[i] );
951 if( ErrorIndex[i] < 2 )
952 ASSERT( acs_abs[p][i] > 0. && acs_sct[p][i] > 0. );
953 if( ErrorIndex[i] < 1 )
956 return lgErrorOccurred;
963 double *normalization)
974 sd->
xx.resize(sd->
nn);
975 sd->
aa.resize(sd->
nn);
976 sd->
rr.resize(sd->
nn);
977 sd->
ww.resize(sd->
nn);
987 for( j=0; j < sd->
nn; j++ )
994 sd->
rr[j] = exp(sd->
rr[j]);
995 sd->
ww[j] *= sd->
rr[j];
1003 *normalization = unity;
1004 sd->
radius *= 1.e-4/unity;
1005 sd->
area *= 4.*
PI*1.e-8/unity;
1006 sd->
vol *= 4./3.*
PI*1.e-12/unity;
1037 const double RATIO_MAX = pow(100.,1./3.);
1044 strcpy( chLine,
" " );
1045 if( strlen(chFile) <= 40 )
1047 strncpy( &chLine[0], chFile, strlen(chFile) );
1051 strncpy( &chLine[0], chFile, 37 );
1052 strncpy( &chLine[37],
"...", 3 );
1055 fprintf(
ioQQQ,
" * >>>> mie_read_opc reading file -- %40s <<<< *\n", chLine );
1058 for(
size_t p=0; p <
gv.ReadRecord.size(); ++p )
1060 if(
gv.ReadRecord[p] == chFile )
1062 fprintf(
ioQQQ,
" File %s has already been read before, was this intended ?\n", chFile );
1066 gv.ReadRecord.push_back( chFile );
1070 nd =
gv.bin.size()-1;
1079 fprintf(
ioQQQ,
" Opacity file %s has obsolete magic number\n",chFile );
1080 fprintf(
ioQQQ,
" I found magic number %ld, but expected %ld on line #%ld\n",
1082 fprintf(
ioQQQ,
" Please recompile this file\n" );
1102 strncpy(
gv.bin[nd]->chDstLab,chLine,(
size_t)
LABELSIZE);
1109 gv.bin[nd]->eec = pow((
double)
gv.bin[nd]->dustp[0],-0.85);
1128 gv.bin[nd]->dustp[4] = 1.;
1133 gv.bin[nd]->eyc = 1./
gv.bin[nd]->AvRadius + 1.e7;
1176 for( nelem=0; nelem <
LIMELM; nelem++ )
1181 gv.bin[nd]->elmAbund[nelem] = RefAbund[nelem];
1184 gv.bin[nd]->AccomCoef[nelem] = 2.*
gv.bin[nd]->atomWeight*
dense.AtomicWeight[nelem]/
1185 pow2(
gv.bin[nd]->atomWeight+
dense.AtomicWeight[nelem]);
1193 lgDefaultQHeat = ( RadiusRatio < RATIO_MAX && !gp.
lgGreyGrain );
1195 gv.bin[nd]->cnv_H_pGR =
gv.bin[nd]->AvVol/
gv.bin[nd]->IntVol;
1196 gv.bin[nd]->cnv_GR_pH = 1./
gv.bin[nd]->cnv_H_pGR;
1204 for( i=0; i < nup; i++ )
1215 double mesh_lo, mesh_hi;
1219 if(
cpu.i().big_endian() )
1220 sscanf( chLine,
"%x %x", &u.i[0], &u.i[1] );
1222 sscanf( chLine,
"%x %x", &u.i[1], &u.i[0] );
1227 if(
cpu.i().big_endian() )
1228 sscanf( chLine,
"%x %x", &u.i[0], &u.i[1] );
1230 sscanf( chLine,
"%x %x", &u.i[1], &u.i[0] );
1233 if( strncmp( md5sum,
continuum.mesh_md5sum.c_str(),
NMD5 ) != 0 ||
1237 fprintf(
ioQQQ,
" Opacity file %s has an incompatible energy grid.\n", chFile );
1238 fprintf(
ioQQQ,
" Please recompile this file using the COMPILE GRAINS command.\n" );
1244 if(
cpu.i().big_endian() )
1245 sscanf( chLine,
"%x %x", &u.i[0], &u.i[1] );
1247 sscanf( chLine,
"%x %x", &u.i[1], &u.i[0] );
1249 gv.bin[nd]->RSFCheck = u.x;
1263 VolTotal =
gv.bin[nd]->IntVol;
1265 for( i=0; i < nbin; i++ )
1273 nd2 =
gv.bin.size()-1;
1276 *
gv.bin[nd2] = *
gv.bin[nd];
1305 gv.bin[nd2]->cnv_H_pGR =
gv.bin[nd2]->AvVol/
gv.bin[nd2]->IntVol;
1306 gv.bin[nd2]->cnv_GR_pH = 1./
gv.bin[nd2]->cnv_H_pGR;
1309 gv.bin[nd2]->Capacity =
1315 gv.bin[nd2]->dustp[4] =
gv.bin[nd2]->IntVol/VolTotal;
1316 for( nelem=0; nelem <
LIMELM; nelem++ )
1317 gv.bin[nd2]->elmAbund[nelem] = RefAbund[nelem]*
gv.bin[nd2]->dustp[4];
1321 for( i=0; i < nbin; i++ )
1327 sprintf(str,
"%02ld",i+1);
1332 for( i=0; i < nbin; i++ )
1335 gv.bin[nd2]->dstab1.resize(nup);
1336 gv.bin[nd2]->pure_sc1.resize(nup);
1337 gv.bin[nd2]->asym.resize(nup);
1338 gv.bin[nd2]->inv_att_len.resize(nup);
1342 for( i=0; i < 5; i++ )
1346 for( i=0; i < nup; i++ )
1349 if( (res = fscanf(io2,
"%le",&anu)) != 1 )
1351 fprintf(
ioQQQ,
" Read failed on %s\n",chFile );
1353 fprintf(
ioQQQ,
" EOF reached prematurely\n" );
1356 for( j=0; j < nbin; j++ )
1359 if( (res = fscanf(io2,
"%le",&
gv.bin[nd2]->dstab1[i])) != 1 )
1361 fprintf(
ioQQQ,
" Read failed on %s\n",chFile );
1363 fprintf(
ioQQQ,
" EOF reached prematurely\n" );
1366 ASSERT(
gv.bin[nd2]->dstab1[i] > 0. );
1371 for( i=0; i < 5; i++ )
1375 for( i=0; i < nup; i++ )
1377 if( (res = fscanf(io2,
"%le",&anu)) != 1 )
1379 fprintf(
ioQQQ,
" Read failed on %s\n",chFile );
1381 fprintf(
ioQQQ,
" EOF reached prematurely\n" );
1384 for( j=0; j < nbin; j++ )
1387 if( (res = fscanf(io2,
"%le",&
gv.bin[nd2]->pure_sc1[i])) != 1 )
1389 fprintf(
ioQQQ,
" Read failed on %s\n",chFile );
1391 fprintf(
ioQQQ,
" EOF reached prematurely\n" );
1394 ASSERT(
gv.bin[nd2]->pure_sc1[i] > 0. );
1399 for( i=0; i < 5; i++ )
1403 for( i=0; i < nup; i++ )
1405 if( (res = fscanf(io2,
"%le",&anu)) != 1 )
1407 fprintf(
ioQQQ,
" Read failed on %s\n",chFile );
1409 fprintf(
ioQQQ,
" EOF reached prematurely\n" );
1412 for( j=0; j < nbin; j++ )
1415 if( (res = fscanf(io2,
"%le",&
gv.bin[nd2]->asym[i])) != 1 )
1417 fprintf(
ioQQQ,
" Read failed on %s\n",chFile );
1419 fprintf(
ioQQQ,
" EOF reached prematurely\n" );
1422 ASSERT(
gv.bin[nd2]->asym[i] > 0. );
1424 gv.bin[nd2]->asym[i] =
min(
gv.bin[nd2]->asym[i],1.);
1429 for( i=0; i < 5; i++ )
1433 for( i=0; i < nup; i++ )
1436 if( (res = fscanf(io2,
"%le %le",&anu,&help)) != 2 )
1438 fprintf(
ioQQQ,
" Read failed on %s\n",chFile );
1440 fprintf(
ioQQQ,
" EOF reached prematurely\n" );
1443 gv.bin[nd]->inv_att_len[i] = (
realnum)help;
1444 ASSERT(
gv.bin[nd]->inv_att_len[i] > 0. );
1446 for( j=1; j < nbin; j++ )
1449 gv.bin[nd2]->inv_att_len[i] =
gv.bin[nd]->inv_att_len[i];
1463 void(*cs_fun)(
double,
const sd_data*,
1492 (*cs_fun)(wavlen,sd,gd,cs_abs,cs_sct,cosb,error);
1511 for( i=0; i < sd->nn; i++ )
1513 sd->cSize = sd->rr[i];
1514 (*cs_fun)(wavlen,sd,gd,&absval,&sctval,&
g,error);
1523 else if( *error == 1 )
1529 *cs_abs += weight*absval;
1530 *cs_sct += weight*sctval;
1531 *cosb += weight*sctval*
g;
1543 *cs_abs /= sd->unity;
1544 *cs_sct /= sd->unity;
1548 fprintf(
ioQQQ,
" insane case for grain size distribution: %d\n" , sd->sdCase );
1554 ASSERT( *cs_abs > 0. && *cs_sct > 0. );
1556 ASSERT( fabs(*cosb) <= 1.+10.*DBL_EPSILON );
1614 ASSERT( frac > 0.-10.*DBL_EPSILON && frac < 1.+10.*DBL_EPSILON );
1615 nre = (1.-frac)*gd->
n[gd->
cAxis][ind].real() + frac*gd->
n[gd->
cAxis][ind+1].real();
1617 nim = (1.-frac)*gd->
n[gd->
cAxis][ind].imag() + frac*gd->
n[gd->
cAxis][ind+1].imag();
1620 ASSERT( fabs(nre-1.-nr1) < 10.*
max(nre,1.)*DBL_EPSILON );
1630 sinpar(nre,nim,x,&qext,&qphase,&qscatt,&ctbrqs,&qback,&iflag);
1638 *cs_abs = area*(qext - qscatt);
1639 *cs_sct = area*qscatt;
1640 *cosb = ctbrqs/qscatt;
1645 if( x >= 100. && sqrt(nr1*nr1+nim*nim) <= 0.001 )
1650 anomal(x,&aqext,&aqabs,&aqphas,&xistar,delta,beta);
1654 *cs_abs = area*aqabs;
1655 *cs_sct = area*(aqext - aqabs);
1669 if( *cs_abs <= 0. || *cs_sct <= 0. )
1671 fprintf(
ioQQQ,
" illegal opacity found: wavl=%.4e micron," , wavlen );
1672 fprintf(
ioQQQ,
" abs_cs=%.2e, sct_cs=%.2e\n" , *cs_abs , *cs_sct );
1673 fprintf(
ioQQQ,
" please check refractive index file...\n" );
1679 if( fabs(*cosb) > 1.+10.*DBL_EPSILON )
1681 fprintf(
ioQQQ,
" illegal asymmetry parameter found: wavl=%.4e micron," , wavlen );
1682 fprintf(
ioQQQ,
" cosb=%.2e\n" , *cosb );
1683 fprintf(
ioQQQ,
" please check refractive index file...\n" );
1711 const double a_xi = 50.e-4;
1712 double xi_PAH, cs_abs1, cs_abs2;
1715 (*pah_fun)(wavl,sd,&gdArr[0],&cs_abs1,cs_sct,cosb,error);
1716 xi_PAH = (1.-q_gra)*
min(1.,
pow3(a_xi/sd->cSize));
1725 mie_cs(wavl,sd,&gdArr[1],&cs_abs2,cs_sct,cosb,error);
1726 *cs_abs = xi_PAH*cs_abs1 + (1.-xi_PAH)*cs_abs2;
1750static const double pah1_strength[7] = { 1.4e-21,1.8e-21,1.2e-20,6.0e-21,4.0e-20,1.9e-20,1.9e-20 };
1751static const double pah1_wlBand[7] = { 3.3, 6.18, 7.7, 8.6, 11.3, 12.0, 13.3 };
1752static const double pah1_width[7] = { 0.024, 0.102, 0.24, 0.168, 0.086, 0.174, 0.174 };
1771 const double p1 = 4.0e-18;
1772 const double p2 = 1.1e-18;
1773 const double p3 = 3.3e-21;
1774 const double p4 = 6.0e-21;
1775 const double p5 = 2.4e-21;
1776 const double wl1a = 5.0;
1777 const double wl1b = 7.0;
1778 const double wl1c = 9.0;
1779 const double wl1d = 9.5;
1780 const double wl2a = 11.0;
1781 const double delwl2 = 0.3;
1783 const double wl2b = wl2a + delwl2;
1784 const double wl2c = 15.0;
1785 const double wl3a = 3.2;
1786 const double wl3b = 3.57;
1787 const double wl3m = (wl3a+wl3b)/2.;
1788 const double wl3sig = 0.1476;
1789 const double x1 = 4.0;
1790 const double x2 = 5.9;
1791 const double x2a =
RYD_INF/1.e4;
1792 const double x3 = 0.1;
1800 double xnh = floor(sqrt(6.*xnc));
1802 double xnhoc = xnh/xnc;
1804 double ftoc3p3 = 100.;
1816 double anu_ev = x/x2a*
EVRYD;
1823 for( j=1; j <= 3; ++j )
1826 csVal1 = (xnh*term1 + xnc*term2)*1.e-18;
1831 cval1 = log(sqrt(xnc)*ftoc3p3/1.2328)/12.2;
1835 term1 = (0.1*x + 0.41)*
pow2(
max(x-
x2,0.));
1859 term2 = exp(cval1*(1. - (
x1/
min(x,
x1))));
1861 term3 = p3*exp(-
pow2(x/x3))*
min(x,x3)/x3;
1863 csVal2 = xnc*((p1*term + p2*term1)*term2 + term3);
1866 if( x2a <= x && x <= 2.*x2a )
1869 double frac =
pow2(2.-x/x2a);
1870 pah1_fun_v = exp((1.-frac)*log(csVal1) + frac*log(csVal2));
1875 pah1_fun_v = csVal1 + csVal2;
1880 if( wl1a <= wavl && wl1d >= wavl )
1883 term = p4*(wavl - wl1a)/(wl1b - wl1a);
1885 term = ( wavl > wl1c ) ? p4*(wl1d - wavl)/(wl1d - wl1c) : p4;
1886 pah1_fun_v += term*xnc;
1888 if( wl2a <= wavl && wl2c >= wavl )
1890 term = ( wavl < wl2b ) ? p5*((wavl - wl2a)/delwl2) : p5*sqrt(1.-
pow2((wavl-wl2a)/(wl2c-wl2a)));
1891 pah1_fun_v += term*xnc;
1893 if( wl3m-10.*wl3sig <= wavl && wavl <= wl3m+10.*wl3sig )
1897 pah1_fun_v += term*xnh;
1901 for( j=0; j < 7; j++ )
1911 if( term1 >= -1. && term1 < -0.5 )
1914 term *= 2. + 2.*term1;
1916 if( term1 >= -0.5 && term1 <= 1.5 )
1918 if( term1 > 1.5 && term1 <= 3. )
1921 term *= 2. - term1*2./3.;
1930 if( term1 >= -2. && term1 < -1. )
1935 if( term1 >= -1. && term1 <= 1. )
1937 if( term1 > 1. && term1 <= 2. )
1943 if( j == 0 || j > 2 )
1945 pah1_fun_v += term*xnc;
1948 *cs_abs = pah1_fun_v;
1951 *cs_sct = 0.1*pah1_fun_v;
1974static const double pah2_wavl[14] = { 0.0722, 0.2175, 3.3, 6.2, 7.7, 8.6, 11.3,
1975 11.9, 12.7, 16.4, 18.3, 21.2, 23.1, 26.0 };
1976static const double pah2_width[14] = { 0.195, 0.217, 0.012, 0.032, 0.091, 0.047, 0.018,
1977 0.025, 0.024, 0.010, 0.036, 0.038, 0.046, 0.69 };
1978static const double pah2n_strength[14] = { 7.97e-13, 1.23e-13, 1.97e-18, 1.96e-19, 6.09e-19, 3.47e-19, 4.27e-18,
1979 7.27e-19, 1.67e-18, 5.52e-20, 6.04e-20, 1.08e-19, 2.78e-20, 1.52e-19 };
1980static const double pah2c_strength[14] = { 7.97e-13, 1.23e-13, 4.47e-19, 1.57e-18, 5.48e-18, 2.42e-18, 4.00e-18,
1981 6.14e-19, 1.49e-18, 5.52e-20, 6.04e-20, 1.08e-19, 2.78e-20, 1.52e-19 };
1997 const double E62 = 3.;
1998 const double E77 = 2.;
1999 const double E86 = 2.;
2010 else if( xnc <= 100. )
2011 xnhoc = 2.5/sqrt(xnc);
2017 double pah2_fun_v = 0.;
2020 double M = ( xnc <= 40. ) ? 0.3*xnc : 0.4*xnc;
2024 cutoff = 1./(3.804/sqrt(
M) + 1.052);
2026 cutoff = 1./(2.282/sqrt(
M) + 0.889);
2027 double y = cutoff/wavl;
2029 double cutoff_fun = atan(1.e3*
pow3(y-1.)/y)/
PI + 0.5;
2031 pah2_fun_v = 34.58*pow( 10., -18.-3.431/x )*cutoff_fun;
2032 for(
int j=2; j < 14; ++j )
2042 strength *= E86*xnhoc;
2043 else if( j == 6 || j == 7 || j == 8 )
2044 strength *= xnhoc/3.;
2052 pah2_fun_v += (1.8687 + 0.1905*x)*1.e-18;
2059 pah2_fun_v += (1.8687 + 0.1905*x +
pow2(y)*(0.4175 + 0.04370*y))*1.e-18;
2064 pah2_fun_v = (((-0.1057*x + 2.950)*x - 24.367)*x + 66.302)*1.e-18;
2070 pah2_fun_v += (-3. + 1.35*x)*1.e-18;
2072 else if( x < 17.26 )
2075 pah2_fun_v = (126.0 - 6.4943*x)*1.e-18;
2085 *cs_abs = pah2_fun_v;
2087 *cs_sct = 0.1*pah2_fun_v;
2110static const double pah3_wavl[30] = { 0.0722, 0.2175, 1.05, 1.26, 1.905, 3.3, 5.27, 5.7, 6.22, 6.69, 7.417, 7.598,
2111 7.85, 8.33, 8.61, 10.68, 11.23, 11.33, 11.99, 12.62, 12.69, 13.48, 14.19,
2112 15.9, 16.45, 17.04, 17.375, 17.87, 18.92, 15. };
2113static const double pah3_width[30] = { 0.195, 0.217, 0.055, 0.11, 0.09, 0.012, 0.034, 0.035, 0.03, 0.07, 0.126,
2114 0.044, 0.053, 0.052, 0.039, 0.02, 0.012, 0.032, 0.045, 0.042, 0.013, 0.04,
2115 0.025, 0.02, 0.014, 0.065, 0.012, 0.016, 0.1, 0.8 };
2116static const double pah3n_strength[30] = { 7.97e-13, 1.23e-13, 0., 0., 0., 3.94e-18, 2.5e-20, 4.e-20, 2.94e-19,
2117 7.35e-20, 2.08e-19, 1.81e-19, 2.19e-19, 6.94e-20, 2.78e-19, 3.e-21,
2118 1.89e-19, 5.2e-19, 2.42e-19, 3.5e-19, 1.3e-20, 8.e-20, 4.5e-21,
2119 4.e-22, 5.e-21, 2.22e-20, 1.1e-21, 6.7e-22, 1.e-21, 5.e-19 };
2120static const double pah3c_strength[30] = { 7.97e-13, 1.23e-13, 2.e-16, 7.8e-17, -1.465e-18, 8.94e-19, 2.e-19,
2121 3.2e-19, 2.35e-18, 5.9e-19, 1.81e-18, 1.63e-18, 1.97e-18, 4.8e-19,
2122 1.94e-18, 3.e-21, 1.77e-19, 4.9e-19, 2.05e-19, 3.1e-19, 1.3e-20,
2123 8.e-20, 4.5e-21, 4.e-22, 5.e-21, 2.22e-20, 1.1e-21, 6.7e-22, 1.7e-21,
2126static const bool pah3_hoc[30] = {
false,
false,
false,
false,
false,
true,
false,
false,
false,
false,
false,
2127 false,
false,
true,
true,
true,
true,
true,
true,
true,
true,
true,
false,
2128 false,
false,
false,
false,
false,
false,
false };
2151 else if( xnc <= 100. )
2152 xnhoc = 2.5/sqrt(xnc);
2159 double pah3_fun_v = ( gd->
charge == 0 ) ? 0. : 3.5*pow(10.,-19.-1.45/x)*exp(-0.1*
pow2(x));
2163 double M = ( xnc <= 40. ) ? 0.3*xnc : 0.4*xnc;
2167 cutoff = 1./(3.804/sqrt(
M) + 1.052);
2169 cutoff = 1./(2.282/sqrt(
M) + 0.889);
2170 double y = cutoff/wavl;
2172 double cutoff_fun = atan(1.e3*
pow3(y-1.)/y)/
PI + 0.5;
2174 pah3_fun_v += 34.58*pow( 10., -18.-3.431/x )*cutoff_fun;
2175 for(
int j=2; j < 30; ++j )
2187 pah3_fun_v += (1.8687 + 0.1905*x)*1.e-18;
2194 pah3_fun_v += (1.8687 + 0.1905*x +
pow2(y)*(0.4175 + 0.04370*y))*1.e-18;
2199 pah3_fun_v += (((-0.1057*x + 2.950)*x - 24.367)*x + 66.302)*1.e-18;
2205 pah3_fun_v += (-3. + 1.35*x)*1.e-18;
2207 else if( x < 17.26 )
2210 pah3_fun_v += (126.0 - 6.4943*x)*1.e-18;
2220 *cs_abs = pah3_fun_v;
2222 *cs_sct = 0.1*pah3_fun_v;
2236 double x = lambda/lambdac;
2238 return 2.e-4/
PI*gamma*lambdac*sigma/(
pow2(x - 1./x) +
pow2(gamma));
2252 double anu =
WAVNRYD/wavl*1.e4;
2263 if( !lgOutOfBounds )
2268 *cs_abs = exp((1.-frac)*log(gd->
opcData[0][ind])+frac*log(gd->
opcData[0][ind+1]));
2271 *cs_sct = exp((1.-frac)*log(gd->
opcData[1][ind])+frac*log(gd->
opcData[1][ind+1]));
2273 *cs_sct = 0.1*(*cs_abs);
2276 a1g = exp((1.-frac)*log(gd->
opcData[2][ind])+frac*log(gd->
opcData[2][ind+1]));
2319 res = pow(size,sd->
a[
ipExp]);
2321 res /= (1. - sd->
a[
ipBeta]*size);
2323 res *= (1. + sd->
a[
ipBeta]*size);
2331 res = exp(-0.5*
pow2(x))/size;
2335 res = exp(-0.5*
pow2(x))/size;
2341 fprintf(
ioQQQ,
" size distribution table has insufficient range\n" );
2342 fprintf(
ioQQQ,
" requested size: %.5f table range %.5f - %.5f\n",
2346 frac = (log(size)-sd->
ln_a[ind])/(sd->
ln_a[ind+1]-sd->
ln_a[ind]);
2347 ASSERT( frac > 0.-10.*DBL_EPSILON && frac < 1.+10.*DBL_EPSILON );
2350 res = exp(res)/
POW4(size);
2354 fprintf(
ioQQQ,
" insane case for grain size distribution: %d\n" , sd->
sdCase );
2384 const double TOLER = 1.e-6;
2389 ASSERT( rel_cutoff > 0. && rel_cutoff < 1. );
2404 f1 = -log(rel_cutoff);
2409 for( i=0; i < 20 && f2 > 0.; ++i )
2422 fprintf(
ioQQQ,
" Could not bracket solution\n" );
2452 vector<double>& invlen,
2453 const char *chString,
2456 bool lgErrorOccurred=
true,
2471 vector<int> ErrorIndex(
rfield.nupper );
2473 for( i=0; i < n; i++ )
2478 lgErrorOccurred =
false;
2481 for( j=0; j < gd->
nAxes; j++ )
2488 lgErrorOccurred =
true;
2493 nim = (1.-frac)*gd->
n[j][ind].imag() + frac*gd->
n[j][ind+1].imag();
2496 InvDep =
PI4*nim/wavlen*1.e4;
2499 invlen[i] += InvDep*gd->
wt[j];
2503 if( lgErrorOccurred )
2505 mie_repair(chString,n,3,3,
rfield.anu,&invlen[0],ErrorIndex,
false,lgWarning);
2523 vector<int>& ErrorIndex,
2543 const long BIG_INTERPOLATION = 10;
2547 lgVerbose = ( chString[0] !=
'\0' );
2549 for( ind1=0; ind1 < n; )
2551 if( ErrorIndex[ind1] == val )
2555 while( ind2 < n && ErrorIndex[ind2] == val )
2559 fprintf(
ioQQQ,
" %s", chString );
2566 lgExtrapolate =
true;
2570 fprintf(
ioQQQ,
" extrapolated below %.4e Ryd\n",anu[i1] );
2573 else if( ind2 == n )
2578 lgExtrapolate =
true;
2582 fprintf(
ioQQQ,
" extrapolated above %.4e Ryd\n",anu[i2] );
2590 lgExtrapolate =
false;
2594 fprintf(
ioQQQ,
" interpolated between %.4e and %.4e Ryd\n",
2597 if( i2-i1-1 > BIG_INTERPOLATION )
2601 fprintf(
ioQQQ,
" ***Warning: extensive interpolation used\n" );
2607 if( i1 < 0 || i2 >= n )
2609 fprintf(
ioQQQ,
" Insufficient data for extrapolation\n" );
2613 xlg1 = log(anu[i1]);
2614 y1lg1 = log(data[i1]);
2617 slp1 =
mie_find_slope(anu,data,ErrorIndex,i1,i2,val,lgVerbose,lgWarning);
2620 xlg2 = log(anu[i2]);
2621 y1lg2 = log(data[i2]);
2622 slp1 = (y1lg2-y1lg1)/(xlg2-xlg1);
2624 if( lgRound && lgExtrapolate && sgn > 0. )
2629 slp1 =
max(slp1,0.);
2632 else if( lgExtrapolate && sgn*slp1 < 0. )
2634 fprintf(
ioQQQ,
" Unphysical value for slope in extrapolation %.6e\n", slp1 );
2635 fprintf(
ioQQQ,
" The most likely cause is that your refractive index or "
2636 "opacity data do not extend to low or high enough frequencies. "
2637 "See Hazy 1 for more details.\n" );
2641 for( j=ind1; j < ind2; j++ )
2643 dx = log(anu[j]) - xlg1;
2644 data[j] = exp(y1lg1 + dx*slp1);
2645 ErrorIndex[j] -= del;
2656 for( j=0; j < n; j++ )
2658 if( ErrorIndex[j] > val-del )
2660 fprintf(
ioQQQ,
" Internal error in mie_repair, index=%ld, val=%d\n",j,ErrorIndex[j] );
2670 const double data[],
2671 const vector<int>& ErrorIndex,
2689 const double LARGE_STDEV = 0.2;
2695 for( i=i1; i <= i2; i++ )
2697 ASSERT( ErrorIndex[i] < val );
2698 ASSERT( anu[i] > 0. && data[i] > 0. );
2710 for( i=i1; i < i2; i++ )
2712 for( j=i+1; j <= i2; j++ )
2714 slp1[k++] = log(data[j]/data[i])/log(anu[j]/anu[i]);
2722 if( slp1[i] > slp1[j] )
2724 double xxx = slp1[i];
2738 s2 +=
pow2(slp1[i]);
2742 stdev = sqrt(stdev);
2745 if( stdev > LARGE_STDEV )
2749 fprintf(
ioQQQ,
" ***Warning: slope for extrapolation may be unreliable\n" );
2761 bool lgLogData=
false;
2794 fprintf(
ioQQQ,
" Refractive index file %s has obsolete magic number\n",chFile );
2796 fprintf(
ioQQQ,
" Please replace this file with an up to date version\n" );
2818 fprintf(
ioQQQ,
" Illegal value for default depletion in %s\n",chFile );
2819 fprintf(
ioQQQ,
" Line #%ld, depl=%14.6e\n",dl,gd->
depl);
2823 for( nelem=0; nelem <
LIMELM; nelem++ )
2836 fprintf(
ioQQQ,
" Illegal value for material type in %s\n",chFile );
2837 fprintf(
ioQQQ,
" Line #%ld, type=%d\n",dl,gd->
matType);
2850 fprintf(
ioQQQ,
" Illegal value for bandgap in %s\n",chFile );
2851 fprintf(
ioQQQ,
" Line #%ld, bandgap=%.4e, work function=%.4e\n",dl,gd->
bandgap,gd->
work);
2852 fprintf(
ioQQQ,
" Bandgap should always be less than work function\n");
2861 fprintf(
ioQQQ,
" Illegal value for thermionic efficiency in %s\n",chFile );
2863 fprintf(
ioQQQ,
" Allowed values are 0. < efficiency <= 1.\n");
2875 if(
nMatch(
"RFI_", chWord ) )
2877 else if(
nMatch(
"OPC_", chWord ) )
2879 else if(
nMatch(
"GREY", chWord ) )
2881 else if(
nMatch(
"PAH1", chWord ) )
2883 else if(
nMatch(
"PH2N", chWord ) )
2885 else if(
nMatch(
"PH2C", chWord ) )
2887 else if(
nMatch(
"PH3N", chWord ) )
2889 else if(
nMatch(
"PH3C", chWord ) )
2893 fprintf(
ioQQQ,
" Illegal keyword in %s\n",chFile );
2894 fprintf(
ioQQQ,
" Line #%ld, value=%s\n",dl,chWord);
2895 fprintf(
ioQQQ,
" Allowed values are: RFI_TBL, OPC_TBL, GREY, PAH1, PH2N, PH2C, PH3N, PH3C\n");
2908 fprintf(
ioQQQ,
" Illegal data code in %s\n",chFile );
2909 fprintf(
ioQQQ,
" Line #%ld, data code=%ld\n",dl,nridf);
2919 fprintf(
ioQQQ,
" Illegal no. of axes in %s\n",chFile );
2920 fprintf(
ioQQQ,
" Line #%ld, number=%ld\n",dl,gd->
nAxes);
2933 if( sscanf( chLine,
"%lf %lf", &gd->
wt[0], &gd->
wt[1] ) != 2 )
2935 fprintf(
ioQQQ,
" Syntax error in %s\n",chFile);
2936 fprintf(
ioQQQ,
" Line #%ld: %s\n",dl,chLine);
2939 if( gd->
wt[0] <= 0. || gd->
wt[1] <= 0. )
2941 fprintf(
ioQQQ,
" Illegal data in %s\n",chFile);
2942 fprintf(
ioQQQ,
" Line #%ld: %s\n",dl,chLine);
2945 total = gd->
wt[0] + gd->
wt[1];
2948 if( sscanf( chLine,
"%lf %lf %lf", &gd->
wt[0], &gd->
wt[1], &gd->
wt[2] ) != 3 )
2950 fprintf(
ioQQQ,
" Syntax error in %s\n",chFile);
2951 fprintf(
ioQQQ,
" Line #%ld: %s\n",dl,chLine);
2954 if( gd->
wt[0] <= 0. || gd->
wt[1] <= 0. || gd->
wt[2] <= 0. )
2956 fprintf(
ioQQQ,
" Illegal data in %s\n",chFile);
2957 fprintf(
ioQQQ,
" Line #%ld: %s\n",dl,chLine);
2960 total = gd->
wt[0] + gd->
wt[1] + gd->
wt[2];
2963 fprintf(
ioQQQ,
" insane case for gd->nAxes: %ld\n", gd->
nAxes );
2967 for( j=0; j < gd->
nAxes; j++ )
2974 if( gd->
ndata[j] < 2 )
2976 fprintf(
ioQQQ,
" Illegal number of data points in %s\n",chFile );
2977 fprintf(
ioQQQ,
" Line #%ld, number=%ld\n",dl,gd->
ndata[j]);
2983 gd->
n[j].resize( gd->
ndata[j] );
2986 for( i=0; i < gd->
ndata[j]; i++ )
2991 if( sscanf( chLine,
"%lf %lf %lf", &gd->
wavlen[j][i], &nr, &ni ) != 3 )
2993 fprintf(
ioQQQ,
" Syntax error in %s\n",chFile);
2994 fprintf(
ioQQQ,
" Line #%ld: %s\n",dl,chLine);
2997 if( gd->
wavlen[j][i] <= 0. )
2999 fprintf(
ioQQQ,
" Illegal value for wavelength in %s\n",chFile);
3000 fprintf(
ioQQQ,
" Line #%ld, wavl=%14.6e\n",dl,gd->
wavlen[j][i]);
3010 fprintf(
ioQQQ,
" Illegal value for wavelength in %s\n",chFile);
3011 fprintf(
ioQQQ,
" Line #%ld, wavl=%14.6e\n",dl,gd->
wavlen[j][i]);
3019 fprintf(
ioQQQ,
" Illegal value for wavelength in %s\n",chFile);
3020 fprintf(
ioQQQ,
" Line #%ld, wavl=%14.6e\n",dl,gd->
wavlen[j][i]);
3032 dftori(&nr,&ni,eps1,eps2);
3033 gd->
nr1[j][i] = nr - 1.;
3040 gd->
nr1[j][i] = nr - 1.;
3043 fprintf(
ioQQQ,
" insane case for nridf: %ld\n", nridf );
3047 gd->
n[j][i] = complex<double>(nr,ni);
3050 if( nr <= 0. || ni < 0. )
3052 fprintf(
ioQQQ,
" Illegal value for refractive index in %s\n",chFile);
3053 fprintf(
ioQQQ,
" Line #%ld, (nr,ni)=(%14.6e,%14.6e)\n",dl,nr,ni);
3056 ASSERT( fabs(nr-1.-gd->
nr1[j][i]) < 10.*nr*DBL_EPSILON );
3071 fprintf(
ioQQQ,
" Illegal no. of data columns in %s\n",chFile );
3072 fprintf(
ioQQQ,
" Line #%ld, number=%ld\n",dl,gd->
nOpcCols);
3080 if(
nMatch(
"LINE", chWord ) )
3082 else if(
nMatch(
"LOG", chWord ) )
3086 fprintf(
ioQQQ,
" Keyword not recognized in %s\n",chFile );
3087 fprintf(
ioQQQ,
" Line #%ld, keyword=%s\n",dl,chWord);
3097 fprintf(
ioQQQ,
" Illegal number of data points in %s\n",chFile );
3098 fprintf(
ioQQQ,
" Line #%ld, number=%ld\n",dl,gd->
nOpcData);
3107 tmp1 = -log10(1.01*DBL_MIN);
3108 tmp2 = log10(0.99*DBL_MAX);
3109 LargestLog =
min(tmp1,tmp2);
3126 if( sscanf( chLine,
"%lf %lf", &gd->
opcAnu[i], &gd->
opcData[0][i] ) != 2 )
3128 fprintf(
ioQQQ,
" Syntax error in %s\n",chFile);
3129 fprintf(
ioQQQ,
" Line #%ld: %s\n",dl,chLine);
3134 if( sscanf( chLine,
"%lf %lf %lf", &gd->
opcAnu[i], &gd->
opcData[0][i],
3137 fprintf(
ioQQQ,
" Syntax error in %s\n",chFile);
3138 fprintf(
ioQQQ,
" Line #%ld: %s\n",dl,chLine);
3143 if( sscanf( chLine,
"%lf %lf %lf %lf", &gd->
opcAnu[i], &gd->
opcData[0][i],
3146 fprintf(
ioQQQ,
" Syntax error in %s\n",chFile);
3147 fprintf(
ioQQQ,
" Line #%ld: %s\n",dl,chLine);
3152 if( sscanf( chLine,
"%lf %lf %lf %lf %lf", &gd->
opcAnu[i], &gd->
opcData[0][i],
3155 fprintf(
ioQQQ,
" Syntax error in %s\n",chFile);
3156 fprintf(
ioQQQ,
" Line #%ld: %s\n",dl,chLine);
3161 fprintf(
ioQQQ,
" insane case for gd->nOpcCols: %ld\n", gd->
nOpcCols );
3168 if( fabs(gd->
opcAnu[i]) > LargestLog )
3170 fprintf(
ioQQQ,
" Illegal value for frequency in %s\n",chFile );
3171 fprintf(
ioQQQ,
" Line #%ld, freq=%14.6e\n",dl,gd->
opcAnu[i] );
3177 if( fabs(gd->
opcData[j][i]) > LargestLog )
3179 fprintf(
ioQQQ,
" Illegal data value in %s\n",chFile );
3180 fprintf(
ioQQQ,
" Line #%ld, value=%14.6e\n",dl,gd->
opcData[j][i] );
3186 if( gd->
opcAnu[i] <= 0. )
3188 fprintf(
ioQQQ,
" Illegal value for frequency in %s\n",chFile );
3189 fprintf(
ioQQQ,
" Line #%ld, freq=%14.6e\n",dl,gd->
opcAnu[i] );
3196 fprintf(
ioQQQ,
" Illegal data value in %s\n",chFile );
3197 fprintf(
ioQQQ,
" Line #%ld, value=%14.6e\n",dl,gd->
opcData[j][i] );
3208 double dataVal = lgLogData ? log10(gd->
opcAnu[1]) : gd->
opcAnu[1];
3209 fprintf(
ioQQQ,
" Illegal value for frequency in %s\n",chFile );
3210 fprintf(
ioQQQ,
" Line #%ld, freq=%14.6e\n",dl,dataVal );
3218 double dataVal = lgLogData ? log10(gd->
opcAnu[i]) : gd->
opcAnu[i];
3219 fprintf(
ioQQQ,
" Illegal value for frequency in %s\n",chFile);
3220 fprintf(
ioQQQ,
" Line #%ld, freq=%14.6e\n",dl,dataVal);
3238 fprintf(
ioQQQ,
" Insane value for gd->rfiType: %d, bailing out\n", gd->
rfiType );
3261 double maxIndex = DBL_MAX,
3269 complex<double> eps_eff(-DBL_MAX,-DBL_MAX);
3286 fprintf(
ioQQQ,
" Mixed grain file %s has obsolete magic number\n",chFile );
3288 fprintf(
ioQQQ,
" Please replace this file with an up to date version\n" );
3297 fprintf(
ioQQQ,
" Illegal value for default depletion in %s\n",chFile );
3298 fprintf(
ioQQQ,
" Line #%ld, depl=%14.6e\n",dl,gd->
depl);
3307 fprintf(
ioQQQ,
" Illegal number of materials in mixed grain file %s\n",chFile );
3308 fprintf(
ioQQQ,
" I found %ld on line #%ld\n",nMaterial,dl );
3309 fprintf(
ioQQQ,
" This number should be at least 2\n" );
3313 vector<double> frac(nMaterial);
3314 vector<double> frac2(nMaterial);
3315 vector<grain_data> gdArr(nMaterial);
3320 for( i=0; i < nMaterial; i++ )
3334 strcpy( chFile2, ++str );
3341 fprintf(
ioQQQ,
" No pair of double quotes was found on line #%ld of file %s\n",dl,chFile );
3342 fprintf(
ioQQQ,
" Please supply the refractive index file name between double quotes\n" );
3349 fprintf(
ioQQQ,
" Input error on line #%ld of file %s\n",dl,chFile );
3350 fprintf(
ioQQQ,
" File %s is not of type RFI_TBL, this is illegal\n",chFile2 );
3355 if( i == nMaterial-1 && gdArr[i].mol_weight == 0. )
3357 fprintf(
ioQQQ,
" Input error on line #%ld of file %s\n",dl,chFile );
3358 fprintf(
ioQQQ,
" Last entry in list of materials is vacuum, this is not allowed\n" );
3359 fprintf(
ioQQQ,
" Please move this entry to an earlier position in the file\n" );
3363 frac2[i] = ( gdArr[i].mol_weight > 0. ) ? frac[i] : 0.;
3366 sumAxes += gdArr[i].nAxes;
3373 if(
nMatch(
"FA00", chWord ) )
3375 else if(
nMatch(
"ST95", chWord ) )
3377 else if(
nMatch(
"BR35", chWord ) )
3381 fprintf(
ioQQQ,
" Keyword not recognized in %s\n",chFile );
3382 fprintf(
ioQQQ,
" Line #%ld, keyword=%s\n",dl,chWord);
3387 for( i=0; i < nMaterial; i++ )
3392 for( nelem=0; nelem <
LIMELM; nelem++ )
3394 gdArr[i].elmAbun[nelem] /= gdArr[i].abun*gdArr[i].depl;
3401 for( nelem=0; nelem <
LIMELM; nelem++ )
3416 for( i=0; i < nMaterial; i++ )
3418 for( k=0; k < gdArr[i].nAxes; k++ )
3420 double wavMin =
min(gdArr[i].wavlen[k][0],gdArr[i].wavlen[k][gdArr[i].ndata[k]-1]);
3421 double wavMax =
max(gdArr[i].wavlen[k][0],gdArr[i].wavlen[k][gdArr[i].ndata[k]-1]);
3422 wavLo =
max(wavLo,wavMin);
3423 wavHi =
min(wavHi,wavMax);
3427 for( nelem=0; nelem <
LIMELM; nelem++ )
3429 gd->
elmAbun[nelem] += frac[i]*gdArr[i].elmAbun[nelem];
3436 gd->
mol_weight += frac[i]*gdArr[i].mol_weight;
3437 gd->
rho += frac[i]*gdArr[i].rho;
3439 if( gdArr[i].mol_weight > 0. )
3446 gd->
work = gdArr[i].work;
3447 gd->
bandgap = gdArr[i].bandgap;
3455 gd->
therm_eff += frac2[i]*gdArr[i].therm_eff;
3458 fprintf(
ioQQQ,
" Insanity in mie_read_mix\n" );
3462 gd->
matType = gdArr[i].matType;
3469 fprintf(
ioQQQ,
" Illegal value for the density: %.3e\n", gd->
rho );
3474 fprintf(
ioQQQ,
" Illegal value for the molecular weight: %.3e\n", gd->
mol_weight );
3477 if( maxIndex <= 0. )
3479 fprintf(
ioQQQ,
" No atoms were found in the grain molecule\n" );
3484 ASSERT( wavLo > 0. && wavHi < DBL_MAX && wavLo < wavHi );
3490 ASSERT( minIndex > 0. && minIndex < DBL_MAX );
3493 wavLo *= 1. + 10.*DBL_EPSILON;
3494 wavHi *= 1. - 10.*DBL_EPSILON;
3498 for( nelem=0; nelem <
LIMELM; nelem++ )
3500 gd->
elmAbun[nelem] /= minIndex;
3504 gd->
abun *= minIndex;
3510 fprintf(
ioQQQ,
"\n The chemical formula of the new grain molecule is: %s\n", chWord );
3511 fprintf(
ioQQQ,
" The abundance wrt H at maximum depletion of this molecule is: %.3e\n",
3513 fprintf(
ioQQQ,
" The abundance wrt H at standard depletion of this molecule is: %.3e\n",
3517 for( nelem=0; nelem <
LIMELM; nelem++ )
3522 vector<double> delta(sumAxes);
3523 vector<double> frdelta(sumAxes);
3524 vector< complex<double> > eps(sumAxes);
3527 for( i=0; i < nMaterial; i++ )
3529 for( k=0; k < gdArr[i].nAxes; k++ )
3531 frdelta[l] = gdArr[i].wt[k]*frac[i];
3532 delta[l] = ( l == 0 ) ? frdelta[l] : delta[l-1] + frdelta[l];
3536 ASSERT( l == sumAxes && fabs(delta[l-1]-1.) < 10.*DBL_EPSILON );
3540 gd->
n[0].resize( gd->
ndata[0] );
3543 wavStep = log(wavHi/wavLo)/(double)(gd->
ndata[0]-1);
3554 for( j=0; j < gd->
ndata[0]; j++ )
3557 complex<double>
a1,
a2,a1c,a2c,a11,a12,a21,a22,ratio;
3559 gd->
wavlen[0][j] = wavLo*exp((
double)j*wavStep);
3563 ratio = eps[0]/eps[1];
3566 a2 = (1.-ratio)*delta[0];
3568 for( l=1; l < sumAxes-1; l++ )
3570 ratio = eps[l]/eps[l+1];
3574 a11 = (ratio+2.)/3.;
3575 a12 = (2.-2.*ratio)/(9.*delta[l]);
3576 a21 = (1.-ratio)*delta[l];
3577 a22 = (2.*ratio+1.)/3.;
3579 a1 = a11*a1c + a12*a2c;
3580 a2 = a21*a1c + a22*a2c;
3587 a21 = eps[sumAxes-1];
3588 a22 = -2./3.*eps[sumAxes-1];
3590 a1 = a11*a1c + a12*a2c;
3591 a2 = a21*a1c + a22*a2c;
3594 dftori(&nre,&nim,ratio.real(),ratio.imag());
3596 gd->
n[0][j] = complex<double>(nre,nim);
3597 gd->
nr1[0][j] = nre-1.;
3602 for( j=0; j < gd->
ndata[0]; j++ )
3604 const double EPS_TOLER = 10.*DBL_EPSILON;
3606 complex<double> eps0;
3608 gd->
wavlen[0][j] = wavLo*exp((
double)j*wavStep);
3617 for( l=0; l < sumAxes; l++ )
3618 eps0 += frdelta[l]*eps[l];
3636 fprintf(
ioQQQ,
" Insanity in mie_read_mix\n" );
3641 dftori(&nre,&nim,eps_eff.real(),eps_eff.imag());
3643 gd->
n[0][j] = complex<double>(nre,nim);
3644 gd->
nr1[0][j] = nre-1.;
3648 fprintf(
ioQQQ,
" Insanity in mie_read_mix\n" );
3659 const vector<grain_data>& gdArr,
3660 vector< complex<double> >& eps)
3668 for( i=0; i < nMaterial; i++ )
3670 for( k=0; k < gdArr[i].nAxes; k++ )
3674 double eps1,eps2,frc,nim,nre;
3676 find_arr(wavlen,gdArr[i].wavlen[k],gdArr[i].ndata[k],&ind,&lgErr);
3678 frc = (wavlen-gdArr[i].wavlen[k][ind])/(gdArr[i].wavlen[k][ind+1]-gdArr[i].wavlen[k][ind]);
3679 ASSERT( frc > 0.-10.*DBL_EPSILON && frc < 1.+10.*DBL_EPSILON );
3680 nre = (1.-frc)*gdArr[i].n[k][ind].real() + frc*gdArr[i].n[k][ind+1].real();
3682 nim = (1.-frc)*gdArr[i].n[k][ind].imag() + frc*gdArr[i].n[k][ind+1].imag();
3684 ritodf(nre,nim,&eps1,&eps2);
3685 eps[l++] = complex<double>(eps1,eps2);
3702 void(*fun)(complex<double>,
const vector<double>&,
const vector< complex<double> >&,
3703 long,complex<double>*,
double*,
double*),
3704 const vector<double>& frdelta,
3705 const vector< complex<double> >& eps,
3713 const double TINY = 1.e-12;
3718 complex<double>
x1,y;
3719 double dudx,dudy,norm2;
3721 (*fun)(
x0,frdelta,eps,sumAxes,&y,&dudx,&dudy);
3725 if( norm2 < TINY*norm(y) )
3727 fprintf(
ioQQQ,
" cnewton - zero divide error\n" );
3731 x1 =
x0 - complex<double>( y.real()*dudx-y.imag()*dudy, y.imag()*dudx+y.real()*dudy )/norm2;
3734 if( fabs(
x0.real()/
x1.real()-1.) + fabs(
x0.imag()/
x1.imag()-1.) < tol )
3742 fprintf(
ioQQQ,
" cnewton did not converge\n" );
3752 const vector<double>& frdelta,
3753 const vector< complex<double> >& eps,
3762 static const double L[4] = { 0., 1./2., 1., 1./3. };
3763 static const double fl[4] = { 5./9., 2./9., 2./9., 1. };
3766 *f = complex<double>(0.,0.);
3769 for( l=0; l < sumAxes; l++ )
3771 complex<double> hlp = eps[l] - x;
3772 double h1 = eps[l].imag()*x.real() - eps[l].real()*x.imag();
3774 for( i=0; i < 4; i++ )
3776 double f1 = fl[i]*frdelta[l];
3777 double xx = ( i < 3 ) ? sin(
PI*frdelta[l]) : cos(
PI*frdelta[l]);
3778 complex<double> f2 = f1*xx*xx;
3779 complex<double> one = x + hlp*L[i];
3780 complex<double> two = f2*hlp/one;
3781 double h2 = norm(one);
3783 *dudx -= f2.real()*(eps[l].real()*
h2 + h1*2.*one.imag()*(1.-L[i]))/
pow2(
h2);
3784 *dudy -= f2.real()*(eps[l].imag()*
h2 - h1*2.*one.real()*(1.-L[i]))/
pow2(
h2);
3794 const vector<double>& frdelta,
3795 const vector< complex<double> >& eps,
3803 static const double L = 1./3.;
3806 *f = complex<double>(0.,0.);
3809 for( l=0; l < sumAxes; l++ )
3811 complex<double> hlp = eps[l] - x;
3812 double h1 = eps[l].imag()*x.real() - eps[l].real()*x.imag();
3813 complex<double> f2 = frdelta[l];
3814 complex<double> one = x + hlp*L;
3815 complex<double> two = f2*hlp/one;
3816 double h2 = norm(one);
3818 *dudx -= f2.real()*(eps[l].real()*
h2 + h1*2.*one.imag()*(1.-L))/
pow2(
h2);
3819 *dudy -= f2.real()*(eps[l].imag()*
h2 - h1*2.*one.real()*(1.-L))/
pow2(
h2);
3829 bool lgTryOverride =
false;
3846 const double FRAC_CUTOFF = 1.e-4;
3847 const double MUL_CO1 = -log(FRAC_CUTOFF);
3848 const double MUL_CO2 = sqrt(MUL_CO1);
3849 const double MUL_CO3 = pow(MUL_CO1,1./3.);
3850 const double MUL_LND = sqrt(-2.*log(FRAC_CUTOFF));
3851 const double MUL_NRM = MUL_LND;
3864 fprintf(
ioQQQ,
" Size distribution file %s has obsolete magic number\n",chFile );
3866 fprintf(
ioQQQ,
" Please replace this file with an up to date version\n" );
3874 if(
nMatch(
"SSIZ", chWord ) )
3878 else if(
nMatch(
"NCAR", chWord ) )
3882 else if(
nMatch(
"POWE", chWord ) )
3886 else if(
nMatch(
"EXP1", chWord ) )
3892 else if(
nMatch(
"EXP2", chWord ) )
3898 else if(
nMatch(
"EXP3", chWord ) )
3904 else if(
nMatch(
"LOGN", chWord ) )
3910 else if(
nMatch(
"NORM", chWord ) )
3915 else if(
nMatch(
"TABL", chWord ) )
3932 fprintf(
ioQQQ,
" Illegal value for grain size\n" );
3933 fprintf(
ioQQQ,
" Grain sizes should be between %.5f and %.0f micron\n",
3935 fprintf(
ioQQQ,
" Line #%ld: %s\n",dl,chLine);
3950 fprintf(
ioQQQ,
" Illegal value for grain size (lower limit)\n" );
3951 fprintf(
ioQQQ,
" Grain sizes should be between %.5f and %.0f micron\n",
3953 fprintf(
ioQQQ,
" Line #%ld: %s\n",dl,chLine);
3963 fprintf(
ioQQQ,
" Illegal value for grain size (upper limit)\n" );
3964 fprintf(
ioQQQ,
" Grain sizes should be between %.5f and %.0f micron\n",
3966 fprintf(
ioQQQ,
" and upper limit should be greater than lower limit\n" );
3967 fprintf(
ioQQQ,
" Line #%ld: %s\n",dl,chLine);
3973 if( sscanf( chLine,
"%lf", &sd->
a[
ipExp] ) != 1 )
3975 fprintf(
ioQQQ,
" Syntax error in %s\n",chFile);
3976 fprintf(
ioQQQ,
" Line #%ld: %s\n",dl,chLine);
4002 if( sscanf( chLine,
"%lf", &sd->
a[
ipExp] ) != 1 )
4004 fprintf(
ioQQQ,
" Syntax error in %s\n",chFile);
4005 fprintf(
ioQQQ,
" Line #%ld: %s\n",dl,chLine);
4011 if( sscanf( chLine,
"%lf", &sd->
a[
ipBeta] ) != 1 )
4013 fprintf(
ioQQQ,
" Syntax error in %s\n",chFile);
4014 fprintf(
ioQQQ,
" Line #%ld: %s\n",dl,chLine);
4027 step_neg = -mul*sd->
a[
ipSLo];
4029 step_pos = mul*sd->
a[
ipSHi];
4030 lgTryOverride =
true;
4045 lgTryOverride =
true;
4058 step_neg = -mul*sd->
a[
ipGSig];
4060 lgTryOverride =
true;
4068 fprintf(
ioQQQ,
" Illegal value for grain size (lower limit)\n" );
4069 fprintf(
ioQQQ,
" Grain sizes should be between %.5f and %.0f micron\n",
4071 fprintf(
ioQQQ,
" Line #%ld: %s\n",dl,chLine);
4081 fprintf(
ioQQQ,
" Illegal value for grain size (upper limit)\n" );
4082 fprintf(
ioQQQ,
" Grain sizes should be between %.5f and %.0f micron\n",
4084 fprintf(
ioQQQ,
" and upper limit should be greater than lower limit\n" );
4085 fprintf(
ioQQQ,
" Line #%ld: %s\n",dl,chLine);
4094 fprintf(
ioQQQ,
" Illegal value for no. of points in table\n" );
4095 fprintf(
ioQQQ,
" Line #%ld: %s\n",dl,chLine);
4104 for( i=0; i < sd->
npts; ++i )
4106 double help1, help2;
4110 if( sscanf( chLine,
"%le %le", &help1, &help2 ) != 2 )
4112 fprintf(
ioQQQ,
" Syntax error in %s\n",chFile);
4113 fprintf(
ioQQQ,
" Line #%ld: %s\n",dl,chLine);
4117 if( help1 <= 0. || help2 <= 0. )
4119 fprintf(
ioQQQ,
" Reading table failed on line #%ld of %s\n",dl,chFile );
4120 fprintf(
ioQQQ,
" Illegal data value %.6e or %.6e\n", help1, help2 );
4124 sd->
ln_a[i] = log(help1);
4127 if( i > 0 && sd->
ln_a[i] <= sd->
ln_a[i-1] )
4129 fprintf(
ioQQQ,
" Reading table failed on line #%ld of %s\n",dl,chFile );
4130 fprintf(
ioQQQ,
" Grain radii should be monotonically increasing\n" );
4140 fprintf(
ioQQQ,
" unimplemented case for grain size distribution in file %s\n" , chFile );
4141 fprintf(
ioQQQ,
" Line #%ld: value %s\n",dl,chWord);
4166 fprintf(
ioQQQ,
" Illegal size limits: lower %.5f and/or upper %.5f\n",
4168 fprintf(
ioQQQ,
" Grain sizes should be between %.5f and %.0f micron\n",
4170 fprintf(
ioQQQ,
" and upper limit should be greater than lower limit\n" );
4171 fprintf(
ioQQQ,
" Please alter the size distribution file\n" );
4182 const char chLine[],
4188 if( sscanf( chLine,
"%ld", data ) != 1 )
4190 fprintf(
ioQQQ,
" Syntax error in %s\n",chFile);
4191 fprintf(
ioQQQ,
" Line #%ld: %s\n",dl,chLine);
4194 if( *data < 0 || (*data == 0 && lgZeroIllegal) )
4196 fprintf(
ioQQQ,
" Illegal data value in %s\n",chFile);
4197 fprintf(
ioQQQ,
" Line #%ld: %ld\n",dl,*data);
4205 const char chLine[],
4212 if( sscanf( chLine,
"%lf", &help ) != 1 )
4214 fprintf(
ioQQQ,
" Syntax error in %s\n",chFile);
4215 fprintf(
ioQQQ,
" Line #%ld: %s\n",dl,chLine);
4219 if( *data < 0. || (*data == 0. && lgZeroIllegal) )
4221 fprintf(
ioQQQ,
" Illegal data value in %s\n",chFile);
4222 fprintf(
ioQQQ,
" Line #%ld: %14.6e\n",dl,*data);
4230 const char chLine[],
4236 if( sscanf( chLine,
"%lf", data ) != 1 )
4238 fprintf(
ioQQQ,
" Syntax error in %s\n",chFile);
4239 fprintf(
ioQQQ,
" Line #%ld: %s\n",dl,chLine);
4242 if( *data < 0. || (*data == 0. && lgZeroIllegal) )
4244 fprintf(
ioQQQ,
" Illegal data value in %s\n",chFile);
4245 fprintf(
ioQQQ,
" Line #%ld: %14.6e\n",dl,*data);
4266 string strWord( chWord );
4267 for( nelem=0; nelem <
LIMELM; nelem++ )
4271 if( chElmName[1] ==
' ' )
4272 chElmName[1] =
'\0';
4273 string::size_type ptr = strWord.find( chElmName );
4274 if( ptr != string::npos )
4276 len = (long)strlen(chElmName);
4278 if( !islower((
unsigned char)chWord[ptr+len]) )
4280 if( isdigit((
unsigned char)chWord[ptr+len]) )
4282 sscanf(chWord+ptr+len,
"%lf",&frac);
4290 elmAbun[nelem] = frac;
4294 *mol_weight += frac*
dense.AtomicWeight[nelem];
4312 for( nelem=0; nelem <
LIMELM; nelem++ )
4314 if( elmAbun[nelem] > 0. )
4317 long index100 =
nint(100.*elmAbun[nelem]);
4320 if( chElmName[1] ==
' ' )
4321 chElmName[1] =
'\0';
4323 if( index100 == 100 )
4324 sprintf( &chWord[len],
"%s", chElmName );
4325 else if( index100%100 == 0 )
4326 sprintf( &chWord[len],
"%s%ld", chElmName, index100/100 );
4329 double xIndex = (double)index100/100.;
4330 sprintf( &chWord[len],
"%s%.2f", chElmName, xIndex );
4332 len = (long)strlen( chWord );
4351 while( chLine[ip] ==
' ' || chLine[ip] ==
'"' )
4354 while( op < n-1 && chLine[ip] !=
' ' && chLine[ip] !=
'"' )
4356 chWord[op++] =
toupper(chLine[ip++]);
4358 chWord[op++] = chLine[ip++];
4379 while( chLine[0] ==
'#' )
4402 fprintf(
ioQQQ,
" Could not read from %s\n",chFile);
4404 fprintf(
ioQQQ,
" EOF reached\n");
4405 fprintf(
ioQQQ,
" This grain data file does not have the expected format.\n");
4428 const vector<double>& x,
4429 const vector<double>& a,
4439 bpa = (xtop+xbot)/2.;
4440 bma = (xtop-xbot)/2.;
4442 for( i=0; i < nn; i++ )
4444 rr[i] = bpa + bma*x[nn-1-i];
4473 const double SAFETY = 5.;
4480 fprintf(
ioQQQ,
" Illegal number of abcissas\n" );
4484 vector<double> c(nn);
4489 for( j=1; j < nn; j++ )
4495 c[j] =
pow2(fj)/((fj-0.5)*(fj+0.5));
4499 for( i=0; i < nn/2; i++ )
4504 xt = 1. - 2.78/(4. +
pow2(fn));
4507 xt = xt - 4.1*(1. + 0.06*(1. - 8./fn))*(1. - xt);
4510 xt = xt - 1.67*(1. + 0.22*(1. - 8./fn))*(x[0] - xt);
4513 xt = 3.*(x[i-1] - x[i-2]) + x[i-3];
4516 for( iter=1; (iter < 20) && (fabs(d) > DBL_EPSILON); iter++ )
4525 for( j=1; j < nn; j++ )
4529 q = 2.*xt*pn - c[j]*pn1;
4530 dq = 2.*xt*dpn - c[j]*dp1 + 2.*pn;
4543 a[i] = cc/(dpn*2.*pn1);
4550 if( fabs(1.-csa) >
SAFETY*fn*DBL_EPSILON )
4552 fprintf(
ioQQQ,
" gauss_legendre failed to converge: delta = %.4e\n", fabs(1.-csa) );
4564 const vector<double>& xa,
4567 bool *lgOutOfBounds)
4580 fprintf(
ioQQQ,
" Invalid array\n");
4586 sgn =
sign3(xa[i3]-xa[i1]);
4589 fprintf(
ioQQQ,
" Ill-ordered array\n");
4593 *lgOutOfBounds = x <
min(xa[0],xa[n-1]) || x >
max(xa[0],xa[n-1]);
4594 if( *lgOutOfBounds )
4601 while( (i3-i1) > 1 )
4603 sgn2 =
sign3(x-xa[i2]);
4687 complex<double> cdum1,
4709 vector< complex<double> > a(
NMXLIM);
4712 ci = complex<double>(0.,1.);
4713 nc = complex<double>(nre,-nim);
4734 xrd = exp(log(x)/3.);
4737 t1 = x + 4.*xrd + 3.;
4741 if( !(x <= 8. || x >= 4200.) )
4742 t1 += 0.05*xrd + 3.;
4753 t4 = x*sqrt(nre*nre+nim*nim);
4776 for( n=0; n < nmx1; n++ )
4779 a[nn-1] = (double)(nn+1)*rrfx - 1./((double)(nn+1)*rrfx+a[nn]);
4784 wn2 = complex<double>(t1,-t2);
4785 wn1 = complex<double>(t2,t1);
4787 tc1 = a[0]*rrf + rx;
4789 sman = (tc1*wn.real() - wn1.real())/(tc1*wn - wn1);
4790 smbn = (tc2*wn.real() - wn1.real())/(tc2*wn - wn1);
4798 nc212 = (nc2-1.)/(nc2+2.);
4802 sman = ci*2.*x3*nc212*(1./3.+x*x*0.2*(nc2-2.)/(nc2+2.)) + 4.*x5*x*nc212*nc212/9.;
4803 smbn = ci*x5*(nc2-1.)/45.;
4812 *qext = 2.*(sman.real() + smbn.real());
4813 *qphase = 2.*(sman.imag() + smbn.imag());
4815 qbck = -2.*(sman - smbn);
4816 *qscat = (norm(sman) + norm(smbn))/.75;
4824 t1 = 2.*(double)n - 1.;
4825 t3 = 2.*(double)n + 1.;
4828 wn = t1*rx*wn1 - wn2;
4831 tc1 = cdum1*rrf + cdum2;
4832 tc2 = cdum1*nc + cdum2;
4833 sman = (tc1*wn.real() - wn1.real())/(tc1*wn - wn1);
4834 smbn = (tc2*wn.real() - wn1.real())/(tc2*wn - wn1);
4838 if( x < xcut && n == 2 )
4841 sman = ci*x5*(nc2-1.)/(15.*(2.*nc2+3.));
4845 eqext = t3*(sman.real() + smbn.real());
4847 eqpha = t3*(sman.imag() + smbn.imag());
4850 eqb = t3*(sman - smbn)*(
double)nsqbk;
4852 tx = norm(sman) + norm(smbn);
4855 t2 = (double)(n - 1);
4858 t2 = (t2*(t5 + 1.))/t5;
4859 ectb = t2*(sman1.real()*sman.real()+sman1.imag()*sman.imag() + smbn1.real()*smbn.real() +
4860 smbn1.imag()*smbn.imag()) +
4861 t4*(sman1.real()*smbn1.real()+sman1.imag()*smbn1.imag());
4869 eqext = fabs(eqext/ *qext);
4870 eqpha = fabs(eqpha/ *qphase);
4871 eqscat = fabs(eqscat/ *qscat);
4872 ectb = ( n == 2 ) ? 0. : fabs(ectb/ *ctbrqs);
4873 eqb = complex<double>( fabs(eqb.real()/qbck.real()), fabs(eqb.imag()/qbck.imag()) );
4875 error =
MAX4(eqext,eqpha,eqscat,ectb);
4877 error1 =
max(eqb.real(),eqb.imag());
4878 if( error < 1.e-07 && error1 < 1.e-04 )
4901 *qback = norm(qbck)*
pow2(rx);
4930 complex<double> cbigk,
4940 *xistar = sqrt(
pow2(xi)+
pow2(xii));
4942 ci = complex<double>(0.,1.);
4943 cw = -complex<double>(xi,xii)*ci;
4945 *qext = 4.*cbigk.real();
4946 *qphase = 4.*cbigk.imag();
4949 *qabs = 2.*cbigk.real();
4956 complex<double> *cbigk)
4964 if( abs(cw) < 1.e-2 )
4969 *cbigk = cw*((1./3.)-cw*((1./8.)-cw*((1./30.)-cw*((1./144.)-cw*((1./840.)-cw*(1./5760.))))));
4973 *cbigk = 0.5 + (exp(-cw)*(1.+cw)-1.)/(cw*cw);
4988 *eps1 = nr*nr - ni*ni;
5004 eps = sqrt(eps2*eps2+eps1*eps1);
5005 *nr = sqrt((eps+eps1)/2.);
5010 *ni = eps2/(2.*(*nr));
const int FILENAME_PATH_LENGTH_2
long nMatch(const char *chKey, const char *chCard)
const char * strstr_s(const char *haystack, const char *needle)
const char * strchr_s(const char *s, int c)
bool fp_equal(sys_float x, sys_float y, int n=3)
char * read_whole_line(char *chLine, int nChar, FILE *ioIN)
NORETURN void TotalInsanity(void)
#define DEBUG_ENTRY(funcname)
double powi(double, long int)
static t_version & Inst()
vector< double > wavlen[NAX]
vector< double > nr1[NAX]
vector< complex< double > > n[NAX]
vector< double > opcData[NDAT]
const multi_geom< d, ALLOC > & clone() const
vector< double > ln_a4dNda
double phfit(long int nz, long int ne, long int is, double e)
void set_version(phfit_version val)
FILE * open_data(const char *fname, const char *mode, access_scheme scheme)
t_elementnames elementnames
static const double TOLER
static const int LABELSUB2
void mie_read_opc(const char *chFile, const GrainPar &gp)
static const long MAGIC_OPC
STATIC void init_eps(double, long, const vector< grain_data > &, vector< complex< double > > &)
static const bool pah3_hoc[30]
STATIC void anomal(double, double *, double *, double *, double *, double, double)
STATIC void dftori(double *, double *, double, double)
STATIC void bigk(complex< double >, complex< double > *)
STATIC void mie_read_long(const char *, const char[], long int *, bool, long int)
STATIC void mie_read_realnum(const char *, const char[], realnum *, bool, long int)
STATIC void mie_integrate(sd_data *, double, double, double *)
void car3_fun(double, const sd_data *, const grain_data[], double *, double *, double *, int *)
static const double pah3_width[30]
static const double pah3n_strength[30]
void mie_write_opc(const char *rfi_file, const char *szd_file, long int nbin)
STATIC void Stognienko(complex< double >, const vector< double > &, const vector< complex< double > > &, long, complex< double > *, double *, double *)
static const long MAGIC_SZD
STATIC void mie_next_line(const char *, FILE *, char *, long int *)
static const int LABELSIZE
static const double LARGEST_GRAIN
STATIC void mie_write_form(const double[], char[])
STATIC complex< double > cnewton(void(*)(complex< double >, const vector< double > &, const vector< complex< double > > &, long, complex< double > *, double *, double *), const vector< double > &, const vector< complex< double > > &, long, complex< double >, double)
static const double pah1_width[7]
STATIC void pah3_fun(double, const sd_data *, const grain_data *, double *, double *, double *, int *)
void gauss_legendre(long int nn, vector< double > &x, vector< double > &a)
static const double pah2n_strength[14]
static const int LABELSUB1
STATIC void tbl_fun(double, const sd_data *, const grain_data *, double *, double *, double *, int *)
STATIC void ritodf(double, double, double *, double *)
STATIC void mie_next_data(const char *, FILE *, char *, long int *)
STATIC void mie_read_word(const char[], char[], long, bool)
static const double pah2_wavl[14]
STATIC void mie_read_szd(const char *, sd_data *)
STATIC void pah2_fun(double, const sd_data *, const grain_data *, double *, double *, double *, int *)
STATIC void mie_repair(const char *, long, int, int, const double[], double[], vector< int > &, bool, bool *)
STATIC double mie_find_slope(const double[], const double[], const vector< int > &, long, long, int, bool, bool *)
double Drude(double, double, double, double)
static const double pah3_wavl[30]
static const double SMALLEST_GRAIN
static const long MAGIC_MIX
STATIC void pah1_fun(double, const sd_data *, const grain_data *, double *, double *, double *, int *)
STATIC void mie_cs_size_distr(double, sd_data *, const grain_data *, void(*)(double, const sd_data *, const grain_data *, double *, double *, double *, int *), double *, double *, double *, int *)
STATIC void mie_read_form(const char *, double[], double *, double *)
static const long MIX_TABLE_SIZE
static const double pah2_width[14]
static const double pah1_strength[7]
void car2_fun(double, const sd_data *, const grain_data[], double *, double *, double *, int *)
STATIC void mie_read_mix(const char *, grain_data *)
STATIC void mie_read_rfi(const char *, grain_data *)
void find_arr(double x, const vector< double > &xa, long int n, long int *ind, bool *lgOutOfBounds)
STATIC double search_limit(double, double, double, sd_data)
STATIC void sinpar(double, double, double, double *, double *, double *, double *, double *, long *)
STATIC void mie_calc_ial(const grain_data *, long, vector< double > &, const char *, bool *)
void car1_fun(double, const sd_data *, const grain_data[], double *, double *, double *, int *)
static const double pah3c_strength[30]
STATIC void ld01_fun(void(*)(double, const sd_data *, const grain_data[], double *, double *, double *, int *), double, double, double, const sd_data *, const grain_data[], double *, double *, double *, int *)
STATIC void mie_auxiliary(sd_data *, const grain_data *, const char *)
static const long MAGIC_RFI
STATIC void mie_read_double(const char *, const char[], double *, bool, long int)
void gauss_init(long int nn, double xbot, double xtop, const vector< double > &x, const vector< double > &a, vector< double > &rr, vector< double > &ww)
STATIC double size_distr(double, const sd_data *)
static const double pah2c_strength[14]
STATIC bool mie_auxiliary2(vector< int > &, multi_arr< double, 2 > &, multi_arr< double, 2 > &, multi_arr< double, 2 > &, long, long)
STATIC void Bruggeman(complex< double >, const vector< double > &, const vector< complex< double > > &, long, complex< double > *, double *, double *)
STATIC void mie_cs(double, const sd_data *, const grain_data *, double *, double *, double *, int *)
static const double pah1_wlBand[7]
static const long LOOP_MAX
static const double SAFETY
double no_atoms(size_t nd)
diatomics h2("h2", 4100., &hmi.H2_total, Yan_H2_CS)
UNUSED const double RYD_INF
UNUSED const double EVRYD
UNUSED const double ELECTRIC_CONST
UNUSED const double ATOMIC_MASS_UNIT
UNUSED const double WAVNRYD
static const unsigned int NMD5