54 static const int NSD = 7;
161 double*,
double*,
int*),
162 double*,
double*,
double*,
int*);
175 long,
long,
int,
bool,
bool*);
180 void(*)(complex<double>,
double[],complex<double>[],
long,complex<double>*,
double*,
double*),
181 double[],complex<double>[],
long,complex<double>,
double);
182 STATIC void Stognienko(complex<double>,
double[],complex<double>[],
long,complex<double>*,
double*,
double*);
183 STATIC void Bruggeman(complex<double>,
double[],complex<double>[],
long,complex<double>*,
double*,
double*);
197 STATIC void sinpar(
double,
double,
double,
double*,
double*,
double*,
198 double*,
double*,
long*);
199 STATIC void anomal(
double,
double*,
double*,
double*,
double*,
double,
double);
200 STATIC void bigk(complex<double>,complex<double>*);
209 const char *szd_file,
248 const long NPTS_TABLE = 100L;
256 for( j=0; j <
NAX; j++ )
264 for( j=0; j <
NDAT; j++ )
269 for( j=0; j <
NAX; j++ )
277 for( j=0; j <
NDAT; j++ )
289 fprintf(
ioQQQ,
" Illegal number of size distribution bins: %ld\n",sd.
nPart );
290 fprintf(
ioQQQ,
" The number should be between 1 and 99.\n" );
295 string rfi_string ( rfi_file );
296 if( rfi_string.find(
".rfi" ) != string::npos )
300 else if( rfi_string.find(
".mix" ) != string::npos )
306 fprintf(
ioQQQ,
" Refractive index file name %s has wrong extention\n", rfi_file );
307 fprintf(
ioQQQ,
" It should have extention .rfi or .mix.\n" );
313 fprintf(
ioQQQ,
" Illegal number of size distribution bins: %ld\n",sd.
nPart );
314 fprintf(
ioQQQ,
" The number should always be 1 for OPC_TABLE files.\n" );
319 fprintf(
ioQQQ,
" Illegal value for the specific density: %.4e\n",gd.
rho );
324 fprintf(
ioQQQ,
" Illegal value for the molecular weight: %.4e\n",gd.
mol_weight );
331 strcpy(chFile,rfi_file);
332 str = strstr(chFile,
".");
337 strcpy(chFile2,szd_file);
338 str = strstr(chFile2,
".");
345 sprintf(ext,
"%02ld",nbin);
346 strcat(strcat(strcat(strcat(strcat(chFile,
"_"),chFile2),
"_"),ext),
".opc");
350 strcat(strcat(strcat(chFile,
"_"),chFile2),
".opc");
353 fprintf(
ioQQQ,
"\n Starting mie_write_opc, output will be written to %s\n\n",chFile );
362 acs_abs = (
double **)
MALLOC(
sizeof(
double *)*(unsigned)sd.
nPart);
363 acs_sct = (
double **)
MALLOC(
sizeof(
double *)*(unsigned)sd.
nPart);
364 a1g = (
double **)
MALLOC(
sizeof(
double *)*(unsigned)sd.
nPart);
367 for( p=0; p < sd.
nPart; p++ )
378 lgErr = lgErr || ( fprintf(fdes,
"# this file was created by Cloudy %s (%s) on %s",
380 lgErr = lgErr || ( fprintf(fdes,
"# ===========================================\n#\n") < 0 );
381 lgErr = lgErr || ( fprintf(fdes,
"%12ld # magic number opacity file\n",
MAGIC_OPC) < 0 );
382 lgErr = lgErr || ( fprintf(fdes,
"%12ld # magic number rfi/mix file\n",gd.
magic) < 0 );
383 lgErr = lgErr || ( fprintf(fdes,
"%12ld # magic number szd file\n",sd.
magic) < 0 );
387 strncpy(hlp1,chFile,(
size_t)(
LABELSUB1+1));
389 str = strstr(hlp1,
"-");
394 strncpy(hlp2,chFile2,(
size_t)(
LABELSUB2+1));
396 str = strstr(hlp2,
"-");
401 strcpy(chGrainLabel,
" ");
406 strcat(strcat(strcat(strcat(chGrainLabel,hlp1),
"-"),hlp2),
"xx");
407 lgErr = lgErr || ( fprintf(fdes,
"%-12.12s # grain type label, xx will be replaced by bin no.\n",
412 strcat(strcat(strcat(chGrainLabel,hlp1),
"-"),hlp2);
413 lgErr = lgErr || ( fprintf(fdes,
"%-12.12s # grain type label\n", chGrainLabel) < 0 );
416 lgErr = lgErr || ( fprintf(fdes,
"%.6e # specific weight (g/cm^3)\n",gd.
rho) < 0 );
417 lgErr = lgErr || ( fprintf(fdes,
"%.6e # molecular weight of grain molecule (amu)\n",gd.
mol_weight) < 0 );
418 lgErr = lgErr || ( fprintf(fdes,
"%.6e # average molecular weight per atom (amu)\n", gd.
atom_weight) < 0 );
419 lgErr = lgErr || ( fprintf(fdes,
"%.6e # abundance of grain molecule at max depletion\n",gd.
abun) < 0 );
420 lgErr = lgErr || ( fprintf(fdes,
"%.6e # depletion efficiency\n",gd.
depl) < 0 );
421 lgErr = lgErr || ( fprintf(fdes,
"%.6e # average grain radius <a^3>/<a^2>, full size distr (cm)\n",
423 lgErr = lgErr || ( fprintf(fdes,
"%.6e # average grain surface area <4pi*a^2>, full size distr (cm^2)\n",
425 lgErr = lgErr || ( fprintf(fdes,
"%.6e # average grain volume <4/3pi*a^3>, full size distr (cm^3)\n",
427 lgErr = lgErr || ( fprintf(fdes,
"%.6e # total grain radius Int(a) per H, full size distr (cm/H)\n",
429 lgErr = lgErr || ( fprintf(fdes,
"%.6e # total grain area Int(4pi*a^2) per H, full size distr (cm^2/H)\n",
431 lgErr = lgErr || ( fprintf(fdes,
"%.6e # total grain vol Int(4/3pi*a^3) per H, full size distr (cm^3/H)\n",
433 lgErr = lgErr || ( fprintf(fdes,
"%.6e # work function (Ryd)\n",gd.
work) < 0 );
434 lgErr = lgErr || ( fprintf(fdes,
"%.6e # gap between valence and conduction band (Ryd)\n",gd.
bandgap) < 0 );
435 lgErr = lgErr || ( fprintf(fdes,
"%.6e # efficiency of thermionic emission\n",gd.
therm_eff) < 0 );
436 lgErr = lgErr || ( fprintf(fdes,
"%.6e # sublimation temperature (K)\n",gd.
subl_temp) < 0 );
437 lgErr = lgErr || ( fprintf(fdes,
"%12d # material type, 1=carbonaceous, 2=silicate\n",gd.
matType) < 0 );
438 lgErr = lgErr || ( fprintf(fdes,
"#\n# abundances of constituent elements rel. to hydrogen\n#\n") < 0 );
440 for( nelem=0; nelem <
LIMELM; nelem++ )
442 lgErr = lgErr || ( fprintf(fdes,
"%.6e # %s\n",gd.
elmAbun[nelem],
448 lgErr = lgErr || ( fprintf(fdes,
"#\n# entire size distribution, amin=%.5f amax=%.5f micron\n",
450 lgErr = lgErr || ( fprintf(fdes,
"#\n%.6e # ratio a_max/a_min in each size bin\n",
452 lgErr = lgErr || ( fprintf(fdes,
"#\n# size distribution function\n#\n") < 0 );
453 lgErr = lgErr || ( fprintf(fdes,
"%12ld # number of table entries\n#\n",NPTS_TABLE+1) < 0 );
454 lgErr = lgErr || ( fprintf(fdes,
"# ============================\n") < 0 );
455 lgErr = lgErr || ( fprintf(fdes,
"# size (micr) a^4*dN/da (cm^3/H)\n#\n") < 0 );
456 for( i=0; i <= NPTS_TABLE; i++ )
462 lgErr = lgErr || ( fprintf(fdes,
"%.6e %.6e\n",radius,a4dNda) < 0 );
467 lgErr = lgErr || ( fprintf(fdes,
"#\n") < 0 );
468 lgErr = lgErr || ( fprintf(fdes,
"%.6e # a_max/a_min = 1 for single sized grain\n", 1. ) < 0 );
469 lgErr = lgErr || ( fprintf(fdes,
"%12ld # no size distribution table\n",0L) < 0 );
472 lgErr = lgErr || ( fprintf(fdes,
"#\n") < 0 );
473 lgErr = lgErr || ( fprintf(fdes,
"%12ld # rfield.nupper\n",
rfield.
nupper) < 0 );
474 lgErr = lgErr || ( fprintf(fdes,
"%12ld # number of size distr. bins\n#\n",sd.
nPart) < 0 );
476 for( p=0; p < sd.
nPart; p++ )
490 volfrac = sd.
vol*frac/volnorm;
491 fprintf(
ioQQQ,
" Starting size bin %ld, amin=%.5f amax=%.5f micron\n",
493 lgErr = lgErr || ( fprintf(fdes,
"# size bin %ld, amin=%.5f amax=%.5f micron\n",
495 lgErr = lgErr || ( fprintf(fdes,
"%.6e # average grain ",3.*sd.
vol/sd.
area) < 0 );
496 lgErr = lgErr || ( fprintf(fdes,
"radius <a^3>/<a^2>, this bin (cm)\n") < 0 );
497 lgErr = lgErr || ( fprintf(fdes,
"%.6e # average ",sd.
area) < 0 );
498 lgErr = lgErr || ( fprintf(fdes,
"grain area <4pi*a^2>, this bin (cm^2)\n") < 0 );
499 lgErr = lgErr || ( fprintf(fdes,
"%.6e # average ",sd.
vol) < 0 );
500 lgErr = lgErr || ( fprintf(fdes,
"grain volume <4/3pi*a^3>, this bin (cm^3)\n") < 0 );
501 lgErr = lgErr || ( fprintf(fdes,
"%.6e # total grain ",sd.
radius*frac/gd.
norm) < 0 );
502 lgErr = lgErr || ( fprintf(fdes,
"radius Int(a) per H, this bin (cm/H)\n") < 0 );
503 lgErr = lgErr || ( fprintf(fdes,
"%.6e # total grain area ",sd.
area*frac/gd.
norm) < 0 );
504 lgErr = lgErr || ( fprintf(fdes,
"Int(4pi*a^2) per H, this bin (cm^2/H)\n") < 0 );
505 lgErr = lgErr || ( fprintf(fdes,
"%.6e # total grain volume ",sd.
vol*frac/gd.
norm) < 0 );
506 lgErr = lgErr || ( fprintf(fdes,
"Int(4/3pi*a^3) per H, this bin (cm^3/H)\n#\n") < 0 );
509 lgErrorOccurred =
false;
527 ErrorIndex[i] =
MAX2(ErrorIndex[i],Error);
528 acs_abs[p][i] += cs_abs*gd.
wt[gd.
cAxis];
529 acs_sct[p][i] += cs_sct*gd.
wt[gd.
cAxis];
530 a1g[p][i] += cs_sct*(1.-cosb)*gd.
wt[gd.
cAxis];
533 if( ErrorIndex[i] > 0 )
535 ErrorIndex[i] =
MIN2(ErrorIndex[i],2);
536 lgErrorOccurred =
true;
539 switch( ErrorIndex[i] )
552 a1g[p][i] /= acs_sct[p][i];
555 fprintf(
ioQQQ,
" Insane value for ErrorIndex: %d\n", ErrorIndex[i] );
561 if( ErrorIndex[i] < 2 )
562 ASSERT( acs_abs[p][i] > 0. && acs_sct[p][i] > 0. );
563 if( ErrorIndex[i] < 1 )
570 ErrorIndex[i] =
MIN2(Error,2);
571 lgErrorOccurred = lgErrorOccurred || ( Error > 0 );
572 acs_abs[p][i] = cs_abs*gd.
norm;
573 acs_sct[p][i] = cs_sct*gd.
norm;
578 acs_abs[p][i] = 1.3121e-23*volfrac*gd.
norm;
579 acs_sct[p][i] = 2.6242e-23*volfrac*gd.
norm;
586 ErrorIndex[i] =
MIN2(Error,2);
587 lgErrorOccurred = lgErrorOccurred || ( Error > 0 );
588 acs_abs[p][i] = cs_abs;
589 acs_sct[p][i] = cs_sct;
593 fprintf(
ioQQQ,
" Insanity in mie_write_opc\n" );
600 if( lgErrorOccurred )
602 strcpy(chString,
"absorption cs");
604 strcpy(chString,
"scattering cs");
606 strcpy(chString,
"asymmetry parameter");
612 acs_abs[p][i] /= gd.
norm;
615 acs_sct[p][i] /= gd.
norm;
620 char name[20] =
"asymxx";
621 sprintf(&name[4],
"%2.2li",p+1);
622 ftmp = fopen(name,
"w");
625 fprintf( ftmp,
"%.6e %.6e\n",
rfield.
anu[i],a1g[p][i]);
634 lgErr = lgErr || ( fprintf(fdes,
"#\n") < 0 );
635 lgErr = lgErr || ( fprintf(fdes,
"# ===========================================\n") < 0 );
636 lgErr = lgErr || ( fprintf(fdes,
"# anu (Ryd) abs_cs_01 (cm^2/H) abs_cs_02.....\n#\n") < 0 );
640 lgErr = lgErr || ( fprintf(fdes,
"%.6e ",
rfield.
anu[i]) < 0 );
641 for( p=0; p < sd.
nPart; p++ )
643 lgErr = lgErr || ( fprintf(fdes,
"%.6e ",acs_abs[p][i]) < 0 );
645 lgErr = lgErr || ( fprintf(fdes,
"\n") < 0 );
648 lgErr = lgErr || ( fprintf(fdes,
"#\n") < 0 );
649 lgErr = lgErr || ( fprintf(fdes,
"# ===========================================\n") < 0 );
650 lgErr = lgErr || ( fprintf(fdes,
"# anu (Ryd) sct_cs_01 (cm^2/H) sct_cs_02.....\n#\n") < 0 );
654 lgErr = lgErr || ( fprintf(fdes,
"%.6e ",
rfield.
anu[i]) < 0 );
655 for( p=0; p < sd.
nPart; p++ )
657 lgErr = lgErr || ( fprintf(fdes,
"%.6e ",acs_sct[p][i]) < 0 );
659 lgErr = lgErr || ( fprintf(fdes,
"\n") < 0 );
662 lgErr = lgErr || ( fprintf(fdes,
"#\n") < 0 );
663 lgErr = lgErr || ( fprintf(fdes,
"# ===========================================\n") < 0 );
664 lgErr = lgErr || ( fprintf(fdes,
"# anu (Ryd) (1-g)_bin_01 (1-g)_bin_02.....\n#\n") < 0 );
668 lgErr = lgErr || ( fprintf(fdes,
"%.6e ",
rfield.
anu[i]) < 0 );
669 for( p=0; p < sd.
nPart; p++ )
671 lgErr = lgErr || ( fprintf(fdes,
"%.6e ",a1g[p][i]) < 0 );
673 lgErr = lgErr || ( fprintf(fdes,
"\n") < 0 );
676 fprintf(
ioQQQ,
" Starting calculation of inverse attenuation length\n" );
677 strcpy(chString,
"inverse attenuation length");
693 fprintf(
ioQQQ,
" mie_write_opc detected unknown ial type: %d\n" , icase );
702 lgErr = lgErr || ( fprintf(fdes,
"#\n") < 0 );
703 lgErr = lgErr || ( fprintf(fdes,
"# ===========================================\n") < 0 );
704 lgErr = lgErr || ( fprintf(fdes,
"# anu (Ryd) inverse attenuation length (cm^-1)\n#\n") < 0 );
708 lgErr = lgErr || ( fprintf(fdes,
"%.6e %.6e\n",
rfield.
anu[i],inv_att_len[i]) < 0 );
715 fprintf(
ioQQQ,
"\n Error writing file: %s\n", chFile );
716 if(
remove(chFile) == 0 )
718 fprintf(
ioQQQ,
" The file has been removed\n" );
724 fprintf(
ioQQQ,
"\n Opacity file %s written succesfully\n\n", chFile );
727 fprintf(
ioQQQ,
"\n !!! Warnings were detected !!!\n\n" );
731 for( p=0; p < sd.
nPart; p++ )
744 if( sd.
ln_a != NULL )
757 for( j=0; j < gd2.
nAxes; j++ )
759 if( gd2.
nr1[j] != NULL )
761 if( gd2.
n[j] != NULL )
763 if( gd2.
wavlen[j] != NULL )
774 for( j=0; j < gd.
nAxes; j++ )
776 if( gd.
nr1[j] != NULL )
778 if( gd.
n[j] != NULL )
780 if( gd.
wavlen[j] != NULL )
799 const double TOLER = 1.e-3;
802 if( strcmp(auxCase,
"init") == 0 )
843 delta = fabs(sd->
vol-oldvol)/sd->
vol;
845 }
while( sd->
nmul <= 1024 && delta > TOLER );
849 fprintf(
ioQQQ,
" could not converge integration of size distribution\n" );
860 fprintf(
ioQQQ,
" insane case for grain size distribution: %d\n" , sd->
sdCase );
865 else if( strcmp(auxCase,
"step") == 0 )
881 step = (amax - amin)/(
double)sd->
nPart;
882 amin = amin + (double)sd->
cPart*step;
883 amax =
MIN2(amax,amin + step);
893 fprintf(
ioQQQ,
" insane case for grain size distribution: %d\n" , sd->
sdCase );
898 else if( strcmp(auxCase,
"cleanup") == 0 )
915 fprintf(
ioQQQ,
" mie_auxiliary called with insane argument: %s\n", auxCase );
925 double *normalization,
937 sd->
xx = (
double *)
MALLOC(
sizeof(
double)*(unsigned)sd->
nn);
938 sd->
aa = (
double *)
MALLOC(
sizeof(
double)*(unsigned)sd->
nn);
939 sd->
rr = (
double *)
MALLOC(
sizeof(
double)*(unsigned)sd->
nn);
940 sd->
ww = (
double *)
MALLOC(
sizeof(
double)*(unsigned)sd->
nn);
950 for( j=0; j < sd->
nn; j++ )
957 sd->
rr[j] = exp(sd->
rr[j]);
958 sd->
ww[j] *= sd->
rr[j];
966 *normalization = unity;
967 sd->
radius *= 1.e-4/unity;
968 sd->
area *= 4.*
PI*1.e-8/unity;
969 sd->
vol *= 4./3.*
PI*1.e-12/unity;
1015 const double RATIO_MAX = pow(100.,1./3.);
1022 strcpy( chLine,
" " );
1023 if( strlen(chFile) <= 40 )
1025 strncpy( &chLine[0], chFile, strlen(chFile) );
1029 strncpy( &chLine[0], chFile, 37 );
1030 strncpy( &chLine[37],
"...", 3 );
1033 fprintf(
ioQQQ,
" * >>>> mie_read_opc reading file -- %40s <<<< *\n", chLine );
1040 fprintf(
ioQQQ,
" File %s has already been read before, was this intended ?\n", chFile );
1061 fprintf(
ioQQQ,
" Opacity file %s has obsolete magic number\n",chFile );
1062 fprintf(
ioQQQ,
" I found magic number %ld, but expected %ld on line #%ld\n",
1064 fprintf(
ioQQQ,
" Please recompile this file\n" );
1155 for( nelem=0; nelem <
LIMELM; nelem++ )
1171 lgDefaultQHeat = ( RadiusRatio < RATIO_MAX && !gp.
lgGreyGrain );
1182 for( i=0; i < nup; i++ )
1208 for( nelem=0; nelem <
LIMELM; nelem++ )
1218 for( i=0; i < nbin; i++ )
1270 for( nelem=0; nelem <
LIMELM; nelem++ )
1280 for( j=0; j < 4; j++ )
1294 for( nelem=0; nelem <
LIMELM; nelem++ )
1300 for( i=0; i < nbin; i++ )
1306 sprintf(str,
"%02ld",i+1);
1311 for( i=0; i < 5; i++ )
1315 for( i=0; i < nup; i++ )
1318 if( (res = fscanf(io2,
"%le",&anu)) != 1 )
1320 fprintf(
ioQQQ,
" Read failed on %s\n",chFile );
1322 fprintf(
ioQQQ,
" EOF reached prematurely\n" );
1325 for( j=0; j < nbin; j++ )
1328 if( (res = fscanf(io2,
"%le",&
gv.
bin[nd2]->
dstab1[i])) != 1 )
1330 fprintf(
ioQQQ,
" Read failed on %s\n",chFile );
1332 fprintf(
ioQQQ,
" EOF reached prematurely\n" );
1340 for( i=0; i < 5; i++ )
1344 for( i=0; i < nup; i++ )
1346 if( (res = fscanf(io2,
"%le",&anu)) != 1 )
1348 fprintf(
ioQQQ,
" Read failed on %s\n",chFile );
1350 fprintf(
ioQQQ,
" EOF reached prematurely\n" );
1353 for( j=0; j < nbin; j++ )
1356 if( (res = fscanf(io2,
"%le",&
gv.
bin[nd2]->
pure_sc1[i])) != 1 )
1358 fprintf(
ioQQQ,
" Read failed on %s\n",chFile );
1360 fprintf(
ioQQQ,
" EOF reached prematurely\n" );
1370 for( i=0; i < 5; i++ )
1374 for( i=0; i < nup; i++ )
1376 if( (res = fscanf(io2,
"%le",&anu)) != 1 )
1378 fprintf(
ioQQQ,
" Read failed on %s\n",chFile );
1380 fprintf(
ioQQQ,
" EOF reached prematurely\n" );
1383 for( j=0; j < nbin; j++ )
1386 if( (res = fscanf(io2,
"%le",&
gv.
bin[nd2]->
asym[i])) != 1 )
1388 fprintf(
ioQQQ,
" Read failed on %s\n",chFile );
1390 fprintf(
ioQQQ,
" EOF reached prematurely\n" );
1398 for( i=0; i < 5; i++ )
1402 for( i=0; i < nup; i++ )
1405 if( (res = fscanf(io2,
"%le %le",&anu,&help)) != 2 )
1407 fprintf(
ioQQQ,
" Read failed on %s\n",chFile );
1409 fprintf(
ioQQQ,
" EOF reached prematurely\n" );
1417 for( j=1; j < nbin; j++ )
1462 (*cs_fun)(wavlen,sd,gd,cs_abs,cs_sct,cosb,error);
1481 for( i=0; i < sd->nn; i++ )
1483 sd->cSize = sd->rr[i];
1484 (*cs_fun)(wavlen,sd,gd,&absval,&sctval,&g,error);
1493 else if( *error == 1 )
1499 *cs_abs += weight*absval;
1500 *cs_sct += weight*sctval;
1501 *cosb += weight*sctval*g;
1513 *cs_abs /= sd->unity;
1514 *cs_sct /= sd->unity;
1518 fprintf(
ioQQQ,
" insane case for grain size distribution: %d\n" , sd->sdCase );
1524 ASSERT( *cs_abs > 0. && *cs_sct > 0. );
1526 ASSERT( fabs(*cosb) <= 1.+10.*DBL_EPSILON );
1584 ASSERT( frac > 0.-10.*DBL_EPSILON && frac < 1.+10.*DBL_EPSILON );
1585 nre = (1.-frac)*gd->
n[gd->
cAxis][ind].real() + frac*gd->
n[gd->
cAxis][ind+1].real();
1587 nim = (1.-frac)*gd->
n[gd->
cAxis][ind].imag() + frac*gd->
n[gd->
cAxis][ind+1].imag();
1590 ASSERT( fabs(nre-1.-nr1) < 10.*
MAX2(nre,1.)*DBL_EPSILON );
1600 sinpar(nre,nim,x,&qext,&qphase,&qscatt,&ctbrqs,&qback,&iflag);
1608 *cs_abs = area*(qext - qscatt);
1609 *cs_sct = area*qscatt;
1610 *cosb = ctbrqs/qscatt;
1615 if( x >= 100. && sqrt(nr1*nr1+nim*nim) <= 0.001 )
1620 anomal(x,&aqext,&aqabs,&aqphas,&xistar,delta,beta);
1624 *cs_abs = area*aqabs;
1625 *cs_sct = area*(aqext - aqabs);
1639 if( *cs_abs <= 0. || *cs_sct <= 0. )
1641 fprintf(
ioQQQ,
" illegal opacity found: wavl=%.4e micron," , wavlen );
1642 fprintf(
ioQQQ,
" abs_cs=%.2e, sct_cs=%.2e\n" , *cs_abs , *cs_sct );
1643 fprintf(
ioQQQ,
" please check refractive index file...\n" );
1649 if( fabs(*cosb) > 1.+10.*DBL_EPSILON )
1651 fprintf(
ioQQQ,
" illegal asymmetry parameter found: wavl=%.4e micron," , wavlen );
1652 fprintf(
ioQQQ,
" cosb=%.2e\n" , *cosb );
1653 fprintf(
ioQQQ,
" please check refractive index file...\n" );
1667 static 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};
1668 static const double pah1_wlBand[7]={3.3, 6.18, 7.7, 8.6, 11.3, 12.0, 13.3};
1669 static const double pah1_width[7]={0.024, 0.102, 0.24, 0.168, 0.086, 0.174, 0.174};
1688 const double p1 = 4.0e-18;
1689 const double p2 = 1.1e-18;
1690 const double p3 = 3.3e-21;
1691 const double p4 = 6.0e-21;
1692 const double p5 = 2.4e-21;
1693 const double wl1a = 5.0;
1694 const double wl1b = 7.0;
1695 const double wl1c = 9.0;
1696 const double wl1d = 9.5;
1697 const double wl2a = 11.0;
1698 const double delwl2 = 0.3;
1700 const double wl2b = wl2a + delwl2;
1701 const double wl2c = 15.0;
1702 const double wl3a = 3.2;
1703 const double wl3b = 3.57;
1704 const double wl3m = (wl3a+wl3b)/2.;
1705 const double wl3sig = 0.1476;
1706 const double x1 = 4.0;
1707 const double x2 = 5.9;
1708 const double x2a =
RYD_INF/1.e4;
1709 const double x3 = 0.1;
1717 double xnh = floor(sqrt(6.*xnc));
1719 double xnhoc = xnh/xnc;
1721 double ftoc3p3 = 100.;
1738 double anu_ev = x/x2a*
EVRYD;
1745 for( j=1; j <= 3; ++j )
1748 csVal1 = (xnh*term1 + xnc*term2)*1.e-18;
1753 cval1 = log(sqrt(xnc)*ftoc3p3/1.2328)/12.2;
1757 term1 = (0.1*x + 0.41)*
POW2(
MAX2(x-x2,0.));
1781 term2 = exp(cval1*(1. - (x1/
MIN2(x,x1))));
1783 term3 = p3*exp(-
POW2(x/x3))*
MIN2(x,x3)/x3;
1785 csVal2 = xnc*((p1*term + p2*term1)*term2 + term3);
1788 if( x2a <= x && x <= 2.*x2a )
1791 double frac =
POW2(2.-x/x2a);
1792 pah1_fun_v = exp((1.-frac)*log(csVal1) + frac*log(csVal2));
1797 pah1_fun_v = csVal1 + csVal2;
1802 if( wl1a <= wavl && wl1d >= wavl )
1805 term = p4*(wavl - wl1a)/(wl1b - wl1a);
1807 term = ( wavl > wl1c ) ? p4*(wl1d - wavl)/(wl1d - wl1c) : p4;
1808 pah1_fun_v += term*xnc;
1810 if( wl2a <= wavl && wl2c >= wavl )
1812 term = ( wavl < wl2b ) ? p5*((wavl - wl2a)/delwl2) : p5*sqrt(1.-
POW2((wavl-wl2a)/(wl2c-wl2a)));
1813 pah1_fun_v += term*xnc;
1815 if( wl3m-10.*wl3sig <= wavl && wavl <= wl3m+10.*wl3sig )
1819 term = 1.1*pah1_strength[0]*exp(-0.5*
POW2((wavl-wl3m)/wl3sig));
1820 pah1_fun_v += term*xnh;
1824 for( j=0; j < 7; j++ )
1826 term1 = (wavl - pah1_wlBand[j])/pah1_width[j];
1834 if( term1 >= -1. && term1 < -0.5 )
1836 term = pah1_strength[j]/(3.*pah1_width[j]);
1837 term *= 2. + 2.*term1;
1839 if( term1 >= -0.5 && term1 <= 1.5 )
1840 term = pah1_strength[j]/(3.*pah1_width[j]);
1841 if( term1 > 1.5 && term1 <= 3. )
1843 term = pah1_strength[j]/(3.*pah1_width[j]);
1844 term *= 2. - term1*2./3.;
1853 if( term1 >= -2. && term1 < -1. )
1855 term = pah1_strength[j]/(3.*pah1_width[j]);
1858 if( term1 >= -1. && term1 <= 1. )
1859 term = pah1_strength[j]/(3.*pah1_width[j]);
1860 if( term1 > 1. && term1 <= 2. )
1862 term = pah1_strength[j]/(3.*pah1_width[j]);
1866 if( j == 0 || j > 2 )
1868 pah1_fun_v += term*xnc;
1871 *cs_abs = pah1_fun_v;
1874 *cs_sct = 0.1*pah1_fun_v;
1892 double anu =
WAVNRYD/wavl*1.e4;
1904 if( !lgOutOfBounds )
1909 *cs_abs = exp((1.-frac)*log(gd->
opcData[0][ind])+frac*log(gd->
opcData[0][ind+1]));
1912 *cs_sct = exp((1.-frac)*log(gd->
opcData[1][ind])+frac*log(gd->
opcData[1][ind+1]));
1914 *cs_sct = 0.1*(*cs_abs);
1917 a1g = exp((1.-frac)*log(gd->
opcData[2][ind])+frac*log(gd->
opcData[2][ind+1]));
1959 res = pow(size,sd->
a[
ipExp]);
1961 res /= (1. - sd->
a[
ipBeta]*size);
1963 res *= (1. + sd->
a[
ipBeta]*size);
1966 if( size > sd->
a[ipBHi] && sd->
a[
ipSHi] > 0. )
1971 res = exp(-0.5*
POW2(x))/size;
1975 res = exp(-0.5*
POW2(x))/size;
1981 fprintf(
ioQQQ,
" size distribution table has insufficient range\n" );
1982 fprintf(
ioQQQ,
" requested size: %.5f table range %.5f - %.5f\n",
1986 frac = (log(size)-sd->
ln_a[ind])/(sd->
ln_a[ind+1]-sd->
ln_a[ind]);
1987 ASSERT( frac > 0.-10.*DBL_EPSILON && frac < 1.+10.*DBL_EPSILON );
1990 res = exp(res)/
POW4(size);
1994 fprintf(
ioQQQ,
" insane case for grain size distribution: %d\n" , sd->
sdCase );
2023 const double TOLER = 1.e-6;
2028 ASSERT( rel_cutoff > 0. && rel_cutoff < 1. );
2043 f1 = -log(rel_cutoff);
2048 for( i=0; i < 20 && f2 > 0.; ++i )
2061 fprintf(
ioQQQ,
" Could not bracket solution\n" );
2066 while( 2.*fabs(x1-x2)/(x1+x2) > TOLER )
2092 const char *chString,
2097 bool lgErrorOccurred=
true,
2115 for( i=0; i < n; i++ )
2120 lgErrorOccurred =
false;
2123 for( j=0; j < gd->
nAxes; j++ )
2130 lgErrorOccurred =
true;
2135 nim = (1.-frac)*gd->
n[j][ind].imag() + frac*gd->
n[j][ind+1].imag();
2138 InvDep =
PI4*nim/wavlen*1.e4;
2141 invlen[i] += InvDep*gd->
wt[j];
2145 if( lgErrorOccurred )
2156 #define NPTS_DERIV 8
2157 #define NPTS_COMB (NPTS_DERIV*(NPTS_DERIV-1)/2)
2186 const long BIG_INTERPOLATION = 10;
2190 lgVerbose = ( chString[0] !=
'\0' );
2192 for( ind1=0; ind1 < n; )
2194 if( ErrorIndex[ind1] == val )
2198 while( ind2 < n && ErrorIndex[ind2] == val )
2202 fprintf(
ioQQQ,
" %s", chString );
2209 lgExtrapolate =
true;
2213 fprintf(
ioQQQ,
" extrapolated below %.4e Ryd\n",anu[i1] );
2216 else if( ind2 == n )
2221 lgExtrapolate =
true;
2225 fprintf(
ioQQQ,
" extrapolated above %.4e Ryd\n",anu[i2] );
2233 lgExtrapolate =
false;
2237 fprintf(
ioQQQ,
" interpolated between %.4e and %.4e Ryd\n",
2240 if( i2-i1-1 > BIG_INTERPOLATION )
2244 fprintf(
ioQQQ,
" ***Warning: extensive interpolation used\n" );
2250 if( i1 < 0 || i2 >= n )
2252 fprintf(
ioQQQ,
" Insufficient data for extrapolation\n" );
2256 xlg1 = log(anu[i1]);
2257 y1lg1 = log(data[i1]);
2260 slp1 =
mie_find_slope(anu,data,ErrorIndex,i1,i2,val,lgVerbose,lgWarning);
2263 xlg2 = log(anu[i2]);
2264 y1lg2 = log(data[i2]);
2265 slp1 = (y1lg2-y1lg1)/(xlg2-xlg1);
2267 if( lgRound && lgExtrapolate && sgn > 0. )
2272 slp1 =
MAX2(slp1,0.);
2275 else if( lgExtrapolate && sgn*slp1 < 0. )
2277 fprintf(
ioQQQ,
" Illegal value for slope in extrapolation %.6e\n", slp1 );
2281 for( j=ind1; j < ind2; j++ )
2283 dx = log(anu[j]) - xlg1;
2284 data[j] = exp(y1lg1 + dx*slp1);
2285 ErrorIndex[j] -= del;
2296 for( j=0; j < n; j++ )
2298 if( ErrorIndex[j] > val-del )
2300 fprintf(
ioQQQ,
" Internal error in mie_repair, index=%ld, val=%d\n",j,ErrorIndex[j] );
2310 const double data[],
2311 const int ErrorIndex[],
2329 const double LARGE_STDEV = 0.2;
2335 for( i=i1; i <= i2; i++ )
2337 ASSERT( ErrorIndex[i] < val );
2338 ASSERT( anu[i] > 0. && data[i] > 0. );
2350 for( i=i1; i < i2; i++ )
2352 for( j=i+1; j <= i2; j++ )
2354 slp1[k++] = log(data[j]/data[i])/log(anu[j]/anu[i]);
2358 for( i=0; i <= NPTS_COMB/2; i++ )
2362 if( slp1[i] > slp1[j] )
2364 double xxx = slp1[i];
2371 slope = ( NPTS_COMB%2 == 1 ) ? slp1[NPTS_COMB/2] : (slp1[NPTS_COMB/2-1] + slp1[NPTS_COMB/2])/2.;
2378 s2 +=
POW2(slp1[i]);
2381 stdev =
MAX2(s2/(
double)NPTS_COMB -
POW2(s1/(
double)NPTS_COMB),0.);
2382 stdev = sqrt(stdev);
2385 for( i=i1; i <= i2; i++ )
2386 printf(
"input: %ld %.4e %.4e\n",i,anu[i],data[i]);
2388 printf(
"%.3f ",slp1[i]);
2390 printf(
"slope %.3f +/- %.3f\n",slope,stdev);
2394 if( stdev > LARGE_STDEV )
2398 fprintf(
ioQQQ,
" ***Warning: slope for extrapolation may be unreliable\n" );
2410 bool lgLogData=
false;
2443 fprintf(
ioQQQ,
" Refractive index file %s has obsolete magic number\n",chFile );
2445 fprintf(
ioQQQ,
" Please replace this file with an up to date version\n" );
2467 fprintf(
ioQQQ,
" Illegal value for default depletion in %s\n",chFile );
2468 fprintf(
ioQQQ,
" Line #%ld, depl=%14.6e\n",dl,gd->
depl);
2472 for( nelem=0; nelem <
LIMELM; nelem++ )
2485 fprintf(
ioQQQ,
" Illegal value for material type in %s\n",chFile );
2486 fprintf(
ioQQQ,
" Line #%ld, type=%d\n",dl,gd->
matType);
2499 fprintf(
ioQQQ,
" Illegal value for bandgap in %s\n",chFile );
2500 fprintf(
ioQQQ,
" Line #%ld, bandgap=%.4e, work function=%.4e\n",dl,gd->
bandgap,gd->
work);
2501 fprintf(
ioQQQ,
" Bandgap should always be less than work function\n");
2510 fprintf(
ioQQQ,
" Illegal value for thermionic efficiency in %s\n",chFile );
2512 fprintf(
ioQQQ,
" Allowed values are 0. < efficiency <= 1.\n");
2524 if(
nMatch(
"RFI_", chWord ) )
2526 else if(
nMatch(
"OPC_", chWord ) )
2528 else if(
nMatch(
"GREY", chWord ) )
2530 else if(
nMatch(
"PAH1", chWord ) )
2534 fprintf(
ioQQQ,
" Illegal keyword in %s\n",chFile );
2535 fprintf(
ioQQQ,
" Line #%ld, value=%s\n",dl,chWord);
2536 fprintf(
ioQQQ,
" Allowed values are: RFI_TBL, OPC_TBL, GREY, PAH1\n");
2549 fprintf(
ioQQQ,
" Illegal data code in %s\n",chFile );
2550 fprintf(
ioQQQ,
" Line #%ld, data code=%ld\n",dl,nridf);
2560 fprintf(
ioQQQ,
" Illegal no. of axes in %s\n",chFile );
2561 fprintf(
ioQQQ,
" Line #%ld, number=%ld\n",dl,gd->
nAxes);
2574 if( sscanf( chLine,
"%lf %lf", &gd->
wt[0], &gd->
wt[1] ) != 2 )
2576 fprintf(
ioQQQ,
" Syntax error in %s\n",chFile);
2577 fprintf(
ioQQQ,
" Line #%ld: %s\n",dl,chLine);
2580 if( gd->
wt[0] <= 0. || gd->
wt[1] <= 0. )
2582 fprintf(
ioQQQ,
" Illegal data in %s\n",chFile);
2583 fprintf(
ioQQQ,
" Line #%ld: %s\n",dl,chLine);
2586 total = gd->
wt[0] + gd->
wt[1];
2589 if( sscanf( chLine,
"%lf %lf %lf", &gd->
wt[0], &gd->
wt[1], &gd->
wt[2] ) != 3 )
2591 fprintf(
ioQQQ,
" Syntax error in %s\n",chFile);
2592 fprintf(
ioQQQ,
" Line #%ld: %s\n",dl,chLine);
2595 if( gd->
wt[0] <= 0. || gd->
wt[1] <= 0. || gd->
wt[2] <= 0. )
2597 fprintf(
ioQQQ,
" Illegal data in %s\n",chFile);
2598 fprintf(
ioQQQ,
" Line #%ld: %s\n",dl,chLine);
2601 total = gd->
wt[0] + gd->
wt[1] + gd->
wt[2];
2604 fprintf(
ioQQQ,
" insane case for gd->nAxes: %ld\n", gd->
nAxes );
2608 for( j=0; j < gd->
nAxes; j++ )
2615 if( gd->
ndata[j] < 2 )
2617 fprintf(
ioQQQ,
" Illegal number of data points in %s\n",chFile );
2618 fprintf(
ioQQQ,
" Line #%ld, number=%ld\n",dl,gd->
ndata[j]);
2625 gd->
n[j] = (complex<double>*)
MALLOC(
sizeof(complex<double>)*(unsigned)gd->
ndata[j]);
2626 gd->
nr1[j] = (
double*)
MALLOC(
sizeof(
double)*(unsigned)gd->
ndata[j]);
2628 for( i=0; i < gd->
ndata[j]; i++ )
2633 if( sscanf( chLine,
"%lf %lf %lf", gd->
wavlen[j]+i, &nr, &ni ) != 3 )
2635 fprintf(
ioQQQ,
" Syntax error in %s\n",chFile);
2636 fprintf(
ioQQQ,
" Line #%ld: %s\n",dl,chLine);
2639 if( gd->
wavlen[j][i] <= 0. )
2641 fprintf(
ioQQQ,
" Illegal value for wavelength in %s\n",chFile);
2642 fprintf(
ioQQQ,
" Line #%ld, wavl=%14.6e\n",dl,gd->
wavlen[j][i]);
2652 fprintf(
ioQQQ,
" Illegal value for wavelength in %s\n",chFile);
2653 fprintf(
ioQQQ,
" Line #%ld, wavl=%14.6e\n",dl,gd->
wavlen[j][i]);
2661 fprintf(
ioQQQ,
" Illegal value for wavelength in %s\n",chFile);
2662 fprintf(
ioQQQ,
" Line #%ld, wavl=%14.6e\n",dl,gd->
wavlen[j][i]);
2674 dftori(&nr,&ni,eps1,eps2);
2675 gd->
nr1[j][i] = nr - 1.;
2682 gd->
nr1[j][i] = nr - 1.;
2685 fprintf(
ioQQQ,
" insane case for nridf: %ld\n", nridf );
2689 gd->
n[j][i] = complex<double>(nr,ni);
2692 if( nr <= 0. || ni < 0. )
2694 fprintf(
ioQQQ,
" Illegal value for refractive index in %s\n",chFile);
2695 fprintf(
ioQQQ,
" Line #%ld, (nr,ni)=(%14.6e,%14.6e)\n",dl,nr,ni);
2698 ASSERT( fabs(nr-1.-gd->
nr1[j][i]) < 10.*nr*DBL_EPSILON );
2713 fprintf(
ioQQQ,
" Illegal no. of data columns in %s\n",chFile );
2714 fprintf(
ioQQQ,
" Line #%ld, number=%ld\n",dl,gd->
nOpcCols);
2722 if(
nMatch(
"LINE", chWord ) )
2724 else if(
nMatch(
"LOG", chWord ) )
2728 fprintf(
ioQQQ,
" Keyword not recognized in %s\n",chFile );
2729 fprintf(
ioQQQ,
" Line #%ld, keyword=%s\n",dl,chWord);
2739 fprintf(
ioQQQ,
" Illegal number of data points in %s\n",chFile );
2740 fprintf(
ioQQQ,
" Line #%ld, number=%ld\n",dl,gd->
nOpcData);
2752 tmp1 = -log10(1.01*DBL_MIN);
2753 tmp2 = log10(0.99*DBL_MAX);
2754 LargestLog =
MIN2(tmp1,tmp2);
2771 if( sscanf( chLine,
"%lf %lf", &gd->
opcAnu[i], &gd->
opcData[0][i] ) != 2 )
2773 fprintf(
ioQQQ,
" Syntax error in %s\n",chFile);
2774 fprintf(
ioQQQ,
" Line #%ld: %s\n",dl,chLine);
2779 if( sscanf( chLine,
"%lf %lf %lf", &gd->
opcAnu[i], &gd->
opcData[0][i],
2782 fprintf(
ioQQQ,
" Syntax error in %s\n",chFile);
2783 fprintf(
ioQQQ,
" Line #%ld: %s\n",dl,chLine);
2788 if( sscanf( chLine,
"%lf %lf %lf %lf", &gd->
opcAnu[i], &gd->
opcData[0][i],
2791 fprintf(
ioQQQ,
" Syntax error in %s\n",chFile);
2792 fprintf(
ioQQQ,
" Line #%ld: %s\n",dl,chLine);
2797 if( sscanf( chLine,
"%lf %lf %lf %lf %lf", &gd->
opcAnu[i], &gd->
opcData[0][i],
2800 fprintf(
ioQQQ,
" Syntax error in %s\n",chFile);
2801 fprintf(
ioQQQ,
" Line #%ld: %s\n",dl,chLine);
2806 fprintf(
ioQQQ,
" insane case for gd->nOpcCols: %ld\n", gd->
nOpcCols );
2813 if( fabs(gd->
opcAnu[i]) > LargestLog )
2815 fprintf(
ioQQQ,
" Illegal value for frequency in %s\n",chFile );
2816 fprintf(
ioQQQ,
" Line #%ld, freq=%14.6e\n",dl,gd->
opcAnu[i] );
2822 if( fabs(gd->
opcData[j][i]) > LargestLog )
2824 fprintf(
ioQQQ,
" Illegal data value in %s\n",chFile );
2825 fprintf(
ioQQQ,
" Line #%ld, value=%14.6e\n",dl,gd->
opcData[j][i] );
2831 if( gd->
opcAnu[i] <= 0. )
2833 fprintf(
ioQQQ,
" Illegal value for frequency in %s\n",chFile );
2834 fprintf(
ioQQQ,
" Line #%ld, freq=%14.6e\n",dl,gd->
opcAnu[i] );
2841 fprintf(
ioQQQ,
" Illegal data value in %s\n",chFile );
2842 fprintf(
ioQQQ,
" Line #%ld, value=%14.6e\n",dl,gd->
opcData[j][i] );
2853 double dataVal = lgLogData ? log10(gd->
opcAnu[1]) : gd->
opcAnu[1];
2854 fprintf(
ioQQQ,
" Illegal value for frequency in %s\n",chFile );
2855 fprintf(
ioQQQ,
" Line #%ld, freq=%14.6e\n",dl,dataVal );
2863 double dataVal = lgLogData ? log10(gd->
opcAnu[i]) : gd->
opcAnu[i];
2864 fprintf(
ioQQQ,
" Illegal value for frequency in %s\n",chFile);
2865 fprintf(
ioQQQ,
" Line #%ld, freq=%14.6e\n",dl,dataVal);
2879 fprintf(
ioQQQ,
" Insane value for gd->rfiType: %d, bailing out\n", gd->
rfiType );
2914 complex<double> *eps,
2915 eps_eff(-DBL_MAX,-DBL_MAX);
2933 fprintf(
ioQQQ,
" Mixed grain file %s has obsolete magic number\n",chFile );
2935 fprintf(
ioQQQ,
" Please replace this file with an up to date version\n" );
2944 fprintf(
ioQQQ,
" Illegal value for default depletion in %s\n",chFile );
2945 fprintf(
ioQQQ,
" Line #%ld, depl=%14.6e\n",dl,gd->
depl);
2954 fprintf(
ioQQQ,
" Illegal number of materials in mixed grain file %s\n",chFile );
2955 fprintf(
ioQQQ,
" I found %ld on line #%ld\n",nMaterial,dl );
2956 fprintf(
ioQQQ,
" This number should be at least 2\n" );
2960 frac = (
double*)
MALLOC(
sizeof(
double)*(unsigned)nMaterial);
2961 frac2 = (
double*)
MALLOC(
sizeof(
double)*(unsigned)nMaterial);
2967 for( i=0; i < nMaterial; i++ )
2972 for( j=0; j <
NAX; j++ )
2974 gdArr[i].
wavlen[j] = NULL;
2975 gdArr[i].
n[j] = NULL;
2976 gdArr[i].
nr1[j] = NULL;
2980 for( j=0; j <
NDAT; j++ )
2992 str = strchr( chLine,
'\"' );
2995 strcpy( chFile2, ++str );
2996 str = strchr( chFile2,
'\"' );
3002 fprintf(
ioQQQ,
" No pair of double quotes was found on line #%ld of file %s\n",dl,chFile );
3003 fprintf(
ioQQQ,
" Please supply the refractive index file name between double quotes\n" );
3010 fprintf(
ioQQQ,
" Input error on line #%ld of file %s\n",dl,chFile );
3011 fprintf(
ioQQQ,
" File %s is not of type RFI_TBL, this is illegal\n",chFile2 );
3016 if( i == nMaterial-1 && gdArr[i].mol_weight == 0. )
3018 fprintf(
ioQQQ,
" Input error on line #%ld of file %s\n",dl,chFile );
3019 fprintf(
ioQQQ,
" Last entry in list of materials is vacuum, this is not allowed\n" );
3020 fprintf(
ioQQQ,
" Please move this entry to an earlier position in the file\n" );
3024 frac2[i] = ( gdArr[i].
mol_weight > 0. ) ? frac[i] : 0.;
3027 sumAxes += gdArr[i].
nAxes;
3034 if(
nMatch(
"FA00", chWord ) )
3036 else if(
nMatch(
"ST95", chWord ) )
3038 else if(
nMatch(
"BR35", chWord ) )
3042 fprintf(
ioQQQ,
" Keyword not recognized in %s\n",chFile );
3043 fprintf(
ioQQQ,
" Line #%ld, keyword=%s\n",dl,chWord);
3048 for( i=0; i < nMaterial; i++ )
3053 for( nelem=0; nelem <
LIMELM; nelem++ )
3062 for( nelem=0; nelem <
LIMELM; nelem++ )
3077 for( i=0; i < nMaterial; i++ )
3079 for( k=0; k < gdArr[i].
nAxes; k++ )
3081 double wavMin =
MIN2(gdArr[i].wavlen[k][0],gdArr[i].wavlen[k][gdArr[i].ndata[k]-1]);
3082 double wavMax =
MAX2(gdArr[i].wavlen[k][0],gdArr[i].wavlen[k][gdArr[i].ndata[k]-1]);
3083 wavLo =
MAX2(wavLo,wavMin);
3084 wavHi =
MIN2(wavHi,wavMax);
3088 for( nelem=0; nelem <
LIMELM; nelem++ )
3098 gd->
rho += frac[i]*gdArr[i].
rho;
3100 if( gdArr[i].mol_weight > 0. )
3119 fprintf(
ioQQQ,
" Insanity in mie_read_mix\n" );
3130 fprintf(
ioQQQ,
" Illegal value for the density: %.3e\n", gd->
rho );
3135 fprintf(
ioQQQ,
" Illegal value for the molecular weight: %.3e\n", gd->
mol_weight );
3138 if( maxIndex <= 0. )
3140 fprintf(
ioQQQ,
" No atoms were found in the grain molecule\n" );
3145 ASSERT( wavLo > 0. && wavHi < DBL_MAX && wavLo < wavHi );
3151 ASSERT( minIndex > 0. && minIndex < DBL_MAX );
3154 wavLo *= 1. + 10.*DBL_EPSILON;
3155 wavHi *= 1. - 10.*DBL_EPSILON;
3159 for( nelem=0; nelem <
LIMELM; nelem++ )
3161 gd->
elmAbun[nelem] /= minIndex;
3165 gd->
abun *= minIndex;
3171 fprintf(
ioQQQ,
"\n The chemical formula of the new grain molecule is: %s\n", chWord );
3172 fprintf(
ioQQQ,
" The abundance wrt H at maximum depletion of this molecule is: %.3e\n",
3174 fprintf(
ioQQQ,
" The abundance wrt H at standard depletion of this molecule is: %.3e\n",
3178 for( nelem=0; nelem <
LIMELM; nelem++ )
3183 delta = (
double*)
MALLOC(
sizeof(
double)*(unsigned)sumAxes);
3184 frdelta = (
double*)
MALLOC(
sizeof(
double)*(unsigned)sumAxes);
3185 eps = (complex<double>*)
MALLOC(
sizeof(complex<double>)*(unsigned)sumAxes);
3188 for( i=0; i < nMaterial; i++ )
3190 for( k=0; k < gdArr[i].
nAxes; k++ )
3192 frdelta[l] = gdArr[i].
wt[k]*frac[i];
3193 delta[l] = ( l == 0 ) ? frdelta[l] : delta[l-1] + frdelta[l];
3197 ASSERT( l == sumAxes && fabs(delta[l-1]-1.) < 10.*DBL_EPSILON );
3202 gd->
n[0] = (complex<double>*)
MALLOC(
sizeof(complex<double>)*(unsigned)gd->
ndata[0]);
3203 gd->
nr1[0] = (
double*)
MALLOC(
sizeof(
double)*(unsigned)gd->
ndata[0]);
3205 wavStep = log(wavHi/wavLo)/(double)(gd->
ndata[0]-1);
3216 for( j=0; j < gd->
ndata[0]; j++ )
3219 complex<double>
a1,
a2,a1c,a2c,a11,a12,a21,a22,ratio;
3221 gd->
wavlen[0][j] = wavLo*exp((
double)j*wavStep);
3225 ratio = eps[0]/eps[1];
3228 a2 = (1.-ratio)*delta[0];
3230 for( l=1; l < sumAxes-1; l++ )
3232 ratio = eps[l]/eps[l+1];
3236 a11 = (ratio+2.)/3.;
3237 a12 = (2.-2.*ratio)/(9.*delta[l]);
3238 a21 = (1.-ratio)*delta[l];
3239 a22 = (2.*ratio+1.)/3.;
3241 a1 = a11*a1c + a12*a2c;
3242 a2 = a21*a1c + a22*a2c;
3249 a21 = eps[sumAxes-1];
3250 a22 = -2./3.*eps[sumAxes-1];
3252 a1 = a11*a1c + a12*a2c;
3253 a2 = a21*a1c + a22*a2c;
3256 dftori(&nre,&nim,ratio.real(),ratio.imag());
3258 gd->
n[0][j] = complex<double>(nre,nim);
3259 gd->
nr1[0][j] = nre-1.;
3264 for( j=0; j < gd->
ndata[0]; j++ )
3266 const double EPS_TOLER = 10.*DBL_EPSILON;
3268 complex<double> eps0;
3270 gd->
wavlen[0][j] = wavLo*exp((
double)j*wavStep);
3279 for( l=0; l < sumAxes; l++ )
3280 eps0 += frdelta[l]*eps[l];
3298 fprintf(
ioQQQ,
" Insanity in mie_read_mix\n" );
3303 dftori(&nre,&nim,eps_eff.real(),eps_eff.imag());
3305 gd->
n[0][j] = complex<double>(nre,nim);
3306 gd->
nr1[0][j] = nre-1.;
3310 fprintf(
ioQQQ,
" Insanity in mie_read_mix\n" );
3316 for( i=0; i < nMaterial; i++ )
3318 for( j=0; j < gdArr[i].
nOpcCols; j++ )
3320 if( gdArr[i].opcData[j] != NULL )
3321 free(gdArr[i].opcData[j]);
3323 if( gdArr[i].opcAnu != NULL )
3324 free(gdArr[i].opcAnu);
3325 for( j=0; j < gdArr[i].
nAxes; j++ )
3327 if( gdArr[i].nr1[j] != NULL )
3328 free(gdArr[i].nr1[j]);
3329 if( gdArr[i].n[j] != NULL )
3330 free(gdArr[i].n[j]);
3331 if( gdArr[i].wavlen[j] != NULL )
3332 free(gdArr[i].wavlen[j]);
3350 complex<double> eps[])
3358 for( i=0; i < nMaterial; i++ )
3360 for( k=0; k < gdArr[i].
nAxes; k++ )
3364 double eps1,eps2,frc,nim,nre;
3366 find_arr(wavlen,gdArr[i].wavlen[k],gdArr[i].ndata[k],&ind,&lgErr);
3368 frc = (wavlen-gdArr[i].
wavlen[k][ind])/(gdArr[i].wavlen[k][ind+1]-gdArr[i].wavlen[k][ind]);
3369 ASSERT( frc > 0.-10.*DBL_EPSILON && frc < 1.+10.*DBL_EPSILON );
3370 nre = (1.-frc)*gdArr[i].n[k][ind].real() + frc*gdArr[i].
n[k][ind+1].real();
3372 nim = (1.-frc)*gdArr[i].n[k][ind].imag() + frc*gdArr[i].
n[k][ind+1].imag();
3374 ritodf(nre,nim,&eps1,&eps2);
3375 eps[l++] = complex<double>(eps1,eps2);
3392 void(*fun)(complex<double>,
double[],complex<double>[],
long,complex<double>*,
double*,
double*),
3394 complex<double> eps[],
3402 const double TINY = 1.e-12;
3407 complex<double>
x1,y;
3408 double dudx,dudy,norm2;
3410 (*fun)(
x0,frdelta,eps,sumAxes,&y,&dudx,&dudy);
3414 if( norm2 < TINY*norm(y) )
3416 fprintf(
ioQQQ,
" cnewton - zero divide error\n" );
3420 x1 = x0 - complex<double>( y.real()*dudx-y.imag()*dudy, y.imag()*dudx+y.real()*dudy )/norm2;
3423 if( fabs(x0.real()/x1.real()-1.) + fabs(x0.imag()/x1.imag()-1.) < tol )
3431 fprintf(
ioQQQ,
" cnewton did not converge\n" );
3442 complex<double> eps[],
3451 static const double L[4] = { 0., 1./2., 1., 1./3. };
3452 static const double fl[4] = { 5./9., 2./9., 2./9., 1. };
3455 *f = complex<double>(0.,0.);
3458 for( l=0; l < sumAxes; l++ )
3460 complex<double> hlp = eps[l] - x;
3461 double h1 = eps[l].imag()*x.real() - eps[l].real()*x.imag();
3463 for( i=0; i < 4; i++ )
3465 double f1 = fl[i]*frdelta[l];
3466 double xx = ( i < 3 ) ? sin(
PI*frdelta[l]) : cos(
PI*frdelta[l]);
3467 complex<double> f2 = f1*xx*xx;
3468 complex<double> one = x + hlp*L[i];
3469 complex<double> two = f2*hlp/one;
3470 double h2 = norm(one);
3472 *dudx -= f2.real()*(eps[l].real()*h2 + h1*2.*one.imag()*(1.-L[i]))/
POW2(h2);
3473 *dudy -= f2.real()*(eps[l].imag()*h2 - h1*2.*one.real()*(1.-L[i]))/
POW2(h2);
3484 complex<double> eps[],
3492 static const double L = 1./3.;
3495 *f = complex<double>(0.,0.);
3498 for( l=0; l < sumAxes; l++ )
3500 complex<double> hlp = eps[l] - x;
3501 double h1 = eps[l].imag()*x.real() - eps[l].real()*x.imag();
3502 complex<double> f2 = frdelta[l];
3503 complex<double> one = x + hlp*L;
3504 complex<double> two = f2*hlp/one;
3505 double h2 = norm(one);
3507 *dudx -= f2.real()*(eps[l].real()*h2 + h1*2.*one.imag()*(1.-L))/
POW2(h2);
3508 *dudy -= f2.real()*(eps[l].imag()*h2 - h1*2.*one.real()*(1.-L))/
POW2(h2);
3518 bool lgTryOverride =
false;
3535 const double FRAC_CUTOFF = 1.e-4;
3536 const double MUL_CO1 = -log(FRAC_CUTOFF);
3537 const double MUL_CO2 = sqrt(MUL_CO1);
3538 const double MUL_CO3 = pow(MUL_CO1,1./3.);
3539 const double MUL_LND = sqrt(-2.*log(FRAC_CUTOFF));
3540 const double MUL_NRM = MUL_LND;
3553 fprintf(
ioQQQ,
" Size distribution file %s has obsolete magic number\n",chFile );
3555 fprintf(
ioQQQ,
" Please replace this file with an up to date version\n" );
3563 if(
nMatch(
"SSIZ", chWord ) )
3567 else if(
nMatch(
"POWE", chWord ) )
3571 else if(
nMatch(
"EXP1", chWord ) )
3577 else if(
nMatch(
"EXP2", chWord ) )
3583 else if(
nMatch(
"EXP3", chWord ) )
3589 else if(
nMatch(
"LOGN", chWord ) )
3595 else if(
nMatch(
"NORM", chWord ) )
3600 else if(
nMatch(
"TABL", chWord ) )
3617 fprintf(
ioQQQ,
" Illegal value for grain size\n" );
3618 fprintf(
ioQQQ,
" Grain sizes should be between %.5f and %.0f micron\n",
3620 fprintf(
ioQQQ,
" Line #%ld: %s\n",dl,chLine);
3630 fprintf(
ioQQQ,
" Illegal value for grain size (lower limit)\n" );
3631 fprintf(
ioQQQ,
" Grain sizes should be between %.5f and %.0f micron\n",
3633 fprintf(
ioQQQ,
" Line #%ld: %s\n",dl,chLine);
3643 fprintf(
ioQQQ,
" Illegal value for grain size (upper limit)\n" );
3644 fprintf(
ioQQQ,
" Grain sizes should be between %.5f and %.0f micron\n",
3646 fprintf(
ioQQQ,
" and upper limit should be greater than lower limit\n" );
3647 fprintf(
ioQQQ,
" Line #%ld: %s\n",dl,chLine);
3653 if( sscanf( chLine,
"%lf", &sd->
a[
ipExp] ) != 1 )
3655 fprintf(
ioQQQ,
" Syntax error in %s\n",chFile);
3656 fprintf(
ioQQQ,
" Line #%ld: %s\n",dl,chLine);
3682 if( sscanf( chLine,
"%lf", &sd->
a[
ipExp] ) != 1 )
3684 fprintf(
ioQQQ,
" Syntax error in %s\n",chFile);
3685 fprintf(
ioQQQ,
" Line #%ld: %s\n",dl,chLine);
3691 if( sscanf( chLine,
"%lf", &sd->
a[
ipBeta] ) != 1 )
3693 fprintf(
ioQQQ,
" Syntax error in %s\n",chFile);
3694 fprintf(
ioQQQ,
" Line #%ld: %s\n",dl,chLine);
3707 step_neg = -mul*sd->
a[
ipSLo];
3709 step_pos = mul*sd->
a[
ipSHi];
3710 lgTryOverride =
true;
3724 step_pos = sd->
a[
ipGCen]*(exp(mul*sd->
a[ipGSig]) - 1.);
3725 lgTryOverride =
true;
3738 step_neg = -mul*sd->
a[
ipGSig];
3740 lgTryOverride =
true;
3748 fprintf(
ioQQQ,
" Illegal value for grain size (lower limit)\n" );
3749 fprintf(
ioQQQ,
" Grain sizes should be between %.5f and %.0f micron\n",
3751 fprintf(
ioQQQ,
" Line #%ld: %s\n",dl,chLine);
3761 fprintf(
ioQQQ,
" Illegal value for grain size (upper limit)\n" );
3762 fprintf(
ioQQQ,
" Grain sizes should be between %.5f and %.0f micron\n",
3764 fprintf(
ioQQQ,
" and upper limit should be greater than lower limit\n" );
3765 fprintf(
ioQQQ,
" Line #%ld: %s\n",dl,chLine);
3774 fprintf(
ioQQQ,
" Illegal value for no. of points in table\n" );
3775 fprintf(
ioQQQ,
" Line #%ld: %s\n",dl,chLine);
3784 for( i=0; i < sd->
npts; ++i )
3786 double help1, help2;
3790 if( sscanf( chLine,
"%le %le", &help1, &help2 ) != 2 )
3792 fprintf(
ioQQQ,
" Syntax error in %s\n",chFile);
3793 fprintf(
ioQQQ,
" Line #%ld: %s\n",dl,chLine);
3797 if( help1 <= 0. || help2 <= 0. )
3799 fprintf(
ioQQQ,
" Reading table failed on line #%ld of %s\n",dl,chFile );
3800 fprintf(
ioQQQ,
" Illegal data value %.6e or %.6e\n", help1, help2 );
3804 sd->
ln_a[i] = log(help1);
3807 if( i > 0 && sd->
ln_a[i] <= sd->
ln_a[i-1] )
3809 fprintf(
ioQQQ,
" Reading table failed on line #%ld of %s\n",dl,chFile );
3810 fprintf(
ioQQQ,
" Grain radii should be monotonically increasing\n" );
3820 fprintf(
ioQQQ,
" unimplemented case for grain size distribution in file %s\n" , chFile );
3821 fprintf(
ioQQQ,
" Line #%ld: value %s\n",dl,chWord);
3846 fprintf(
ioQQQ,
" Illegal size limits: lower %.5f and/or upper %.5f\n",
3848 fprintf(
ioQQQ,
" Grain sizes should be between %.5f and %.0f micron\n",
3850 fprintf(
ioQQQ,
" and upper limit should be greater than lower limit\n" );
3851 fprintf(
ioQQQ,
" Please alter the size distribution file\n" );
3862 const char chLine[],
3868 if( sscanf( chLine,
"%ld", data ) != 1 )
3870 fprintf(
ioQQQ,
" Syntax error in %s\n",chFile);
3871 fprintf(
ioQQQ,
" Line #%ld: %s\n",dl,chLine);
3874 if( *data < 0 || (*data == 0 && lgZeroIllegal) )
3876 fprintf(
ioQQQ,
" Illegal data value in %s\n",chFile);
3877 fprintf(
ioQQQ,
" Line #%ld: %ld\n",dl,*data);
3885 const char chLine[],
3892 if( sscanf( chLine,
"%lf", &help ) != 1 )
3894 fprintf(
ioQQQ,
" Syntax error in %s\n",chFile);
3895 fprintf(
ioQQQ,
" Line #%ld: %s\n",dl,chLine);
3899 if( *data < 0. || (*data == 0. && lgZeroIllegal) )
3901 fprintf(
ioQQQ,
" Illegal data value in %s\n",chFile);
3902 fprintf(
ioQQQ,
" Line #%ld: %14.6e\n",dl,*data);
3910 const char chLine[],
3916 if( sscanf( chLine,
"%lf", data ) != 1 )
3918 fprintf(
ioQQQ,
" Syntax error in %s\n",chFile);
3919 fprintf(
ioQQQ,
" Line #%ld: %s\n",dl,chLine);
3922 if( *data < 0. || (*data == 0. && lgZeroIllegal) )
3924 fprintf(
ioQQQ,
" Illegal data value in %s\n",chFile);
3925 fprintf(
ioQQQ,
" Line #%ld: %14.6e\n",dl,*data);
3946 string strWord( chWord );
3947 for( nelem=0; nelem <
LIMELM; nelem++ )
3951 if( chElmName[1] ==
' ' )
3952 chElmName[1] =
'\0';
3953 string::size_type ptr = strWord.find( chElmName );
3954 if( ptr != string::npos )
3956 len = (long)strlen(chElmName);
3958 if( !islower((
unsigned char)chWord[ptr+len]) )
3960 if( isdigit((
unsigned char)chWord[ptr+len]) )
3962 sscanf(chWord+ptr+len,
"%lf",&frac);
3970 elmAbun[nelem] = frac;
3977 if( *no_atoms == 0. )
3992 for( nelem=0; nelem <
LIMELM; nelem++ )
3994 if( elmAbun[nelem] > 0. )
3997 long index100 =
nint(100.*elmAbun[nelem]);
4000 if( chElmName[1] ==
' ' )
4001 chElmName[1] =
'\0';
4003 if( index100 == 100 )
4004 sprintf( &chWord[len],
"%s", chElmName );
4005 else if( index100%100 == 0 )
4006 sprintf( &chWord[len],
"%s%ld", chElmName, index100/100 );
4009 double xIndex = (double)index100/100.;
4010 sprintf( &chWord[len],
"%s%.2f", chElmName, xIndex );
4012 len = (long)strlen( chWord );
4031 while( chLine[ip] ==
' ' || chLine[ip] ==
'"' )
4034 while( op < n-1 && chLine[ip] !=
' ' && chLine[ip] !=
'"' )
4036 chWord[op++] = (char)toupper(chLine[ip++]);
4038 chWord[op++] = chLine[ip++];
4059 while( chLine[0] ==
'#' )
4065 str = strstr(chLine,
"#");
4082 fprintf(
ioQQQ,
" Could not read from %s\n",chFile);
4084 fprintf(
ioQQQ,
" EOF reached\n");
4118 bpa = (xtop+xbot)/2.;
4119 bma = (xtop-xbot)/2.;
4121 for( i=0; i < nn; i++ )
4123 rr[i] = bpa + bma*x[nn-1-i];
4153 const double SAFETY = 5.;
4160 fprintf(
ioQQQ,
" Illegal number of abcissas\n" );
4164 c = (
double *)
MALLOC(
sizeof(
double)*(unsigned)nn);
4169 for( j=1; j < nn; j++ )
4175 c[j] =
POW2(fj)/((fj-0.5)*(fj+0.5));
4179 for( i=0; i < nn/2; i++ )
4184 xt = 1. - 2.78/(4. +
POW2(fn));
4187 xt = xt - 4.1*(1. + 0.06*(1. - 8./fn))*(1. - xt);
4190 xt = xt - 1.67*(1. + 0.22*(1. - 8./fn))*(x[0] - xt);
4193 xt = 3.*(x[i-1] - x[i-2]) + x[i-3];
4196 for( iter=1; (iter < 20) && (fabs(d) > DBL_EPSILON); iter++ )
4205 for( j=1; j < nn; j++ )
4209 q = 2.*xt*pn - c[j]*pn1;
4210 dq = 2.*xt*dpn - c[j]*dp1 + 2.*pn;
4223 a[i] = cc/(dpn*2.*pn1);
4230 if( fabs(1.-csa) > SAFETY*fn*DBL_EPSILON )
4232 fprintf(
ioQQQ,
" gauss_legendre failed to converge: delta = %.4e\n",fabs(1.-csa) );
4248 bool *lgOutOfBounds)
4261 fprintf(
ioQQQ,
" Invalid array\n");
4267 sgn =
sign3(xa[i3]-xa[i1]);
4270 fprintf(
ioQQQ,
" Ill-ordered array\n");
4274 *lgOutOfBounds = x <
MIN2(xa[0],xa[n-1]) || x >
MAX2(xa[0],xa[n-1]);
4275 if( *lgOutOfBounds )
4282 while( (i3-i1) > 1 )
4284 sgn2 =
sign3(x-xa[i2]);
4334 #define NMXLIM 16000
4372 complex<double> cdum1,
4394 a = (complex<double>*)
MALLOC(
sizeof(complex<double>)*(unsigned)
NMXLIM);
4397 ci = complex<double>(0.,1.);
4398 nc = complex<double>(nre,-nim);
4419 xrd = exp(log(x)/3.);
4422 t1 = x + 4.*xrd + 3.;
4426 if( !(x <= 8. || x >= 4200.) )
4427 t1 += 0.05*xrd + 3.;
4438 t4 = x*sqrt(nre*nre+nim*nim);
4461 for( n=0; n < nmx1; n++ )
4464 a[nn-1] = (double)(nn+1)*rrfx - 1./((double)(nn+1)*rrfx+a[nn]);
4469 wn2 = complex<double>(t1,-t2);
4470 wn1 = complex<double>(t2,t1);
4472 tc1 = a[0]*rrf + rx;
4474 sman = (tc1*wn.real() - wn1.real())/(tc1*wn - wn1);
4475 smbn = (tc2*wn.real() - wn1.real())/(tc2*wn - wn1);
4483 nc212 = (nc2-1.)/(nc2+2.);
4487 sman = ci*2.*x3*nc212*(1./3.+x*x*0.2*(nc2-2.)/(nc2+2.)) + 4.*x5*x*nc212*nc212/9.;
4488 smbn = ci*x5*(nc2-1.)/45.;
4497 *qext = 2.*(sman.real() + smbn.real());
4498 *qphase = 2.*(sman.imag() + smbn.imag());
4500 qbck = -2.*(sman - smbn);
4501 *qscat = (norm(sman) + norm(smbn))/.75;
4509 t1 = 2.*(double)n - 1.;
4510 t3 = 2.*(double)n + 1.;
4513 wn = t1*rx*wn1 - wn2;
4516 tc1 = cdum1*rrf + cdum2;
4517 tc2 = cdum1*nc + cdum2;
4518 sman = (tc1*wn.real() - wn1.real())/(tc1*wn - wn1);
4519 smbn = (tc2*wn.real() - wn1.real())/(tc2*wn - wn1);
4523 if( x < xcut && n == 2 )
4526 sman = ci*x5*(nc2-1.)/(15.*(2.*nc2+3.));
4530 eqext = t3*(sman.real() + smbn.real());
4532 eqpha = t3*(sman.imag() + smbn.imag());
4535 eqb = t3*(sman - smbn)*(
double)nsqbk;
4537 tx = norm(sman) + norm(smbn);
4540 t2 = (double)(n - 1);
4543 t2 = (t2*(t5 + 1.))/t5;
4544 ectb = t2*(sman1.real()*sman.real()+sman1.imag()*sman.imag() + smbn1.real()*smbn.real() +
4545 smbn1.imag()*smbn.imag()) +
4546 t4*(sman1.real()*smbn1.real()+sman1.imag()*smbn1.imag());
4554 eqext = fabs(eqext/ *qext);
4555 eqpha = fabs(eqpha/ *qphase);
4556 eqscat = fabs(eqscat/ *qscat);
4557 ectb = ( n == 2 ) ? 0. : fabs(ectb/ *ctbrqs);
4558 eqb = complex<double>( fabs(eqb.real()/qbck.real()), fabs(eqb.imag()/qbck.imag()) );
4560 error =
MAX4(eqext,eqpha,eqscat,ectb);
4562 error1 =
MAX2(eqb.real(),eqb.imag());
4563 if( error < 1.e-07 && error1 < 1.e-04 )
4586 *qback = norm(qbck)*
POW2(rx);
4617 complex<double> cbigk,
4627 *xistar = sqrt(
POW2(xi)+
POW2(xii));
4629 ci = complex<double>(0.,1.);
4630 cw = -complex<double>(xi,xii)*ci;
4632 *qext = 4.*cbigk.real();
4633 *qphase = 4.*cbigk.imag();
4636 *qabs = 2.*cbigk.real();
4643 complex<double> *cbigk)
4651 if( abs(cw) < 1.e-2 )
4656 *cbigk = cw*((1./3.)-cw*((1./8.)-cw*((1./30.)-cw*((1./144.)-cw*((1./840.)-cw*(1./5760.))))));
4660 *cbigk = 0.5 + (exp(-cw)*(1.+cw)-1.)/(cw*cw);
4675 *eps1 = nr*nr - ni*ni;
4691 eps = sqrt(eps2*eps2+eps1*eps1);
4692 *nr = sqrt((eps+eps1)/2.);
4697 *ni = eps2/(2.*(*nr));