00001
00002
00003 #include "cddefines.h"
00004 #include "continuum.h"
00005 #include "thirdparty.h"
00006 #include "elementnames.h"
00007 #include "physconst.h"
00008 #include "dense.h"
00009 #include "called.h"
00010 #include "version.h"
00011 #include "grainvar.h"
00012 #include "rfield.h"
00013 #include "atmdat.h"
00014 #include "grains.h"
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038 static const long MAGIC_RFI = 1030103L;
00039 static const long MAGIC_SZD = 2010403L;
00040 static const long MAGIC_OPC = 3100827L;
00041 static const long MAGIC_MIX = 4030103L;
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051 static const double SMALLEST_GRAIN = 0.0001*(1.-10.*DBL_EPSILON);
00052 static const double LARGEST_GRAIN = 10.*(1.+10.*DBL_EPSILON);
00053
00054
00055 static const int NSD = 7;
00056
00057
00058
00059 static const int ipSize = 0;
00060 static const int ipBLo = 0;
00061 static const int ipBHi = 1;
00062 static const int ipExp = 2;
00063 static const int ipBeta = 3;
00064 static const int ipSLo = 4;
00065 static const int ipSHi = 5;
00066 static const int ipAlpha = 6;
00067 static const int ipGCen = 2;
00068 static const int ipGSig = 3;
00070
00071 typedef enum {
00072 RFI_TABLE, OPC_TABLE, OPC_GREY, OPC_PAH1, OPC_PAH2N, OPC_PAH2C
00073 } rfi_type;
00074
00075
00076 typedef enum {
00077 FARAFONOV00, STOGNIENKO95, BRUGGEMAN35
00078 } emt_type;
00079
00080
00081 typedef enum {
00082 SD_ILLEGAL, SD_SINGLE_SIZE, SD_POWERLAW, SD_EXP_CUTOFF1, SD_EXP_CUTOFF2,
00083 SD_EXP_CUTOFF3, SD_LOG_NORMAL, SD_LIN_NORMAL, SD_TABLE, SD_NR_CARBON
00084 } sd_type;
00085
00086 class sd_data {
00087 void p_clear1()
00088 {
00089 xx.clear();
00090 aa.clear();
00091 rr.clear();
00092 ww.clear();
00093 ln_a.clear();
00094 ln_a4dNda.clear();
00095 }
00096 public:
00097 double a[NSD];
00098 double lim[2];
00099 double clim[2];
00100 vector<double> xx;
00101 vector<double> aa;
00102 vector<double> rr;
00103 vector<double> ww;
00104 double unity;
00105 double unity_bin;
00106 double cSize;
00107 double radius;
00108 double area;
00109 double vol;
00110 vector<double> ln_a;
00111 vector<double> ln_a4dNda;
00112 sd_type sdCase;
00113 long int nCarbon;
00114 long int magic;
00115 long int cPart;
00116 long int nPart;
00117 long int nmul;
00118 long int nn;
00119 long int npts;
00120 bool lgLogScale;
00121 void clear()
00122 {
00123 p_clear1();
00124 }
00125 ~sd_data()
00126 {
00127 p_clear1();
00128 }
00129 };
00130
00131
00132 static const int NAX = 3;
00133 static const int NDAT = 4;
00134
00135 class grain_data {
00136 void p_clear0()
00137 {
00138 nAxes = 0;
00139 nOpcCols = 0;
00140 }
00141 void p_clear1()
00142 {
00143 for( int j=0; j < NAX; j++ )
00144 {
00145 wavlen[j].clear();
00146 n[j].clear();
00147 nr1[j].clear();
00148 }
00149 opcAnu.clear();
00150 for( int j=0; j < NDAT; j++ )
00151 opcData[j].clear();
00152 }
00153 public:
00154 vector<double> wavlen[NAX];
00155 vector< complex<double> > n[NAX];
00156 vector<double> nr1[NAX];
00157 vector<double> opcAnu;
00158 vector<double> opcData[NDAT];
00159 double wt[NAX];
00160 double abun;
00161 double depl;
00162 double elmAbun[LIMELM];
00163 double mol_weight;
00164 double atom_weight;
00165 double rho;
00166 double norm;
00167 double work;
00168 double bandgap;
00169 double therm_eff;
00170 double subl_temp;
00171 long int magic;
00172 long int cAxis;
00173 long int nAxes;
00174 long int ndata[NAX];
00175 long int nOpcCols;
00176 long int nOpcData;
00177 long int charge;
00178 mat_type matType;
00179 rfi_type rfiType;
00180 void clear()
00181 {
00182 p_clear1();
00183 p_clear0();
00184 }
00185 grain_data()
00186 {
00187 p_clear0();
00188 }
00189 ~grain_data()
00190 {
00191 p_clear1();
00192 }
00193 };
00194
00195 static const int WORDLEN = 5;
00196
00197
00198 static const int LABELSUB1 = 3;
00199 static const int LABELSUB2 = 5;
00200 static const int LABELSIZE = LABELSUB1 + LABELSUB2 + 4;
00201
00202
00203 static const long MIX_TABLE_SIZE = 2000L;
00204
00205 STATIC void mie_auxiliary(sd_data*,const grain_data*,const char*);
00206 STATIC bool mie_auxiliary2(vector<int>&,multi_arr<double,2>&,
00207 multi_arr<double,2>&,multi_arr<double,2>&,long,long);
00208 STATIC void mie_integrate(sd_data*,double,double,double*);
00209 STATIC void mie_cs_size_distr(double,sd_data*,const grain_data*,
00210 void(*)(double,const sd_data*,const grain_data*,
00211 double*,double*,double*,int*),
00212 double*,double*,double*,int*);
00213 STATIC void mie_cs(double,const sd_data*,const grain_data*,double*,double*,
00214 double*,int*);
00215 STATIC void car1_fun(double,const sd_data*,const grain_data[],double*,double*,
00216 double*,int*);
00217 STATIC void pah1_fun(double,const sd_data*,const grain_data*,double*,double*,
00218 double*,int*);
00219 STATIC void ld01_fun(double,const sd_data*,const grain_data[],double*,double*,
00220 double*,int*);
00221 STATIC void pah2_fun(double,const sd_data*,const grain_data*,double*,double*,
00222 double*,int*);
00223 inline double Drude(double,double,double,double);
00224 STATIC void tbl_fun(double,const sd_data*,const grain_data*,double*,double*,
00225 double*,int*);
00226 STATIC double size_distr(double,const sd_data*);
00227 STATIC double search_limit(double,double,double,sd_data);
00228 STATIC void mie_calc_ial(const grain_data*,long,vector<double>&,const char*,bool*);
00229 STATIC void mie_repair(const char*,long,int,int,const realnum[],double[],vector<int>&,
00230 bool,bool*);
00231 STATIC double mie_find_slope(const realnum[],const double[],const vector<int>&,
00232 long,long,int,bool,bool*);
00233 STATIC void mie_read_rfi(const char*,grain_data*);
00234 STATIC void mie_read_mix(const char*,grain_data*);
00235 STATIC void init_eps(double,long,const vector<grain_data>&,vector< complex<double> >&);
00236 STATIC complex<double> cnewton(
00237 void(*)(complex<double>,const vector<double>&,const vector< complex<double> >&,
00238 long,complex<double>*,double*,double*),
00239 const vector<double>&,const vector< complex<double> >&,long,complex<double>,double);
00240 STATIC void Stognienko(complex<double>,const vector<double>&,const vector< complex<double> >&,
00241 long,complex<double>*,double*,double*);
00242 STATIC void Bruggeman(complex<double>,const vector<double>&,const vector< complex<double> >&,
00243 long,complex<double>*,double*,double*);
00244 STATIC void mie_read_szd(const char*,sd_data*);
00245 STATIC void mie_read_long(const char*,const char[],long int*,bool,long int);
00246 STATIC void mie_read_realnum(const char*,const char[],realnum*,bool,long int);
00247 STATIC void mie_read_double(const char*,const char[],double*,bool,long int);
00248 STATIC void mie_read_form(const char*,double[],double*,double*);
00249 STATIC void mie_write_form(const double[],char[]);
00250 STATIC void mie_read_word(const char[],char[],long,bool);
00251 STATIC void mie_next_data(const char*,FILE*,char*,long int*);
00252 STATIC void mie_next_line(const char*,FILE*,char*,long int*);
00253
00254
00255
00256
00257 STATIC void sinpar(double,double,double,double*,double*,double*,
00258 double*,double*,long*);
00259 STATIC void anomal(double,double*,double*,double*,double*,double,double);
00260 STATIC void bigk(complex<double>,complex<double>*);
00261 STATIC void ritodf(double,double,double*,double*);
00262 STATIC void dftori(double*,double*,double,double);
00263
00264
00265 void mie_write_opc( const char *rfi_file,
00266 const char *szd_file,
00267 long int nbin)
00268 {
00269 int Error = 0;
00270 bool lgErr,
00271 lgErrorOccurred,
00272 lgWarning;
00273 long int i,
00274 nelem,
00275 p;
00276 double cosb,
00277 cs_abs,
00278 cs_sct,
00279 volfrac,
00280 volnorm,
00281 wavlen;
00282 char chGrainLabel[LABELSIZE+1],
00283 ext[3],
00284 chFile[FILENAME_PATH_LENGTH_2],
00285 chFile2[FILENAME_PATH_LENGTH_2],
00286 hlp1[LABELSUB1+2],
00287 hlp2[LABELSUB2+2],
00288 *str,
00289 chString[FILENAME_PATH_LENGTH_2];
00290 time_t timer;
00291 FILE *fdes;
00292
00293
00294
00295 const long NPTS_TABLE = 100L;
00296
00297 DEBUG_ENTRY( "mie_write_opc()" );
00298
00299 sd_data sd;
00300
00301 mie_read_szd( szd_file , &sd );
00302
00303 sd.nPart = ( sd.sdCase == SD_SINGLE_SIZE || sd.sdCase == SD_NR_CARBON ) ? 1 : nbin;
00304 if( sd.nPart <= 0 || sd.nPart >= 100 )
00305 {
00306 fprintf( ioQQQ, " Illegal number of size distribution bins: %ld\n", sd.nPart );
00307 fprintf( ioQQQ, " The number should be between 1 and 99.\n" );
00308 cdEXIT(EXIT_FAILURE);
00309 }
00310 sd.lgLogScale = true;
00311
00312 vector<grain_data> gdArr(2);
00313 grain_data& gd = gdArr[0];
00314 grain_data& gd2 = gdArr[1];
00315
00316 string rfi_string ( rfi_file );
00317 if( rfi_string.find( ".rfi" ) != string::npos )
00318 {
00319 mie_read_rfi( rfi_file , &gd );
00320 }
00321 else if( rfi_string.find( ".mix" ) != string::npos )
00322 {
00323 mie_read_mix( rfi_file , &gd );
00324 }
00325 else
00326 {
00327 fprintf( ioQQQ, " Refractive index file name %s has wrong extention\n", rfi_file );
00328 fprintf( ioQQQ, " It should have extention .rfi or .mix.\n" );
00329 cdEXIT(EXIT_FAILURE);
00330 }
00331
00332 if( gd.rfiType == OPC_TABLE && sd.nPart > 1 )
00333 {
00334 fprintf( ioQQQ, " Illegal number of size distribution bins: %ld\n", sd.nPart );
00335 fprintf( ioQQQ, " The number should always be 1 for OPC_TABLE files.\n" );
00336 cdEXIT(EXIT_FAILURE);
00337 }
00338 if( gd.rho <= 0. )
00339 {
00340 fprintf( ioQQQ, " Illegal value for the specific density: %.4e\n", gd.rho );
00341 cdEXIT(EXIT_FAILURE);
00342 }
00343 if( gd.mol_weight <= 0. )
00344 {
00345 fprintf( ioQQQ, " Illegal value for the molecular weight: %.4e\n", gd.mol_weight );
00346 cdEXIT(EXIT_FAILURE);
00347 }
00348
00349 lgWarning = false;
00350
00351
00352 strcpy(chFile,rfi_file);
00353 str = strstr_s(chFile,".");
00354
00355 if( str != NULL )
00356 *str = '\0';
00357
00358 strcpy(chFile2,szd_file);
00359 str = strstr_s(chFile2,".");
00360
00361 if( str != NULL )
00362 *str = '\0';
00363
00364 if( sd.sdCase != SD_SINGLE_SIZE && sd.sdCase != SD_NR_CARBON )
00365 {
00366 sprintf(ext,"%02ld",nbin);
00367 strcat(strcat(strcat(strcat(strcat(chFile,"_"),chFile2),"_"),ext),".opc");
00368 }
00369 else
00370 {
00371 strcat(strcat(strcat(chFile,"_"),chFile2),".opc");
00372 }
00373
00374 mie_auxiliary(&sd,&gd,"init");
00375
00376
00377 gd.norm = sd.vol*gd.rho/(ATOMIC_MASS_UNIT*gd.mol_weight*gd.abun*gd.depl);
00378 volnorm = sd.vol;
00379 volfrac = 1.;
00380
00381 multi_arr<double,2> acs_abs( sd.nPart, rfield.nupper );
00382 multi_arr<double,2> acs_sct( acs_abs.clone() );
00383 multi_arr<double,2> a1g( acs_abs.clone() );
00384 vector<double> inv_att_len( rfield.nupper );
00385
00386 fprintf( ioQQQ, "\n Starting mie_write_opc, output will be written to %s\n\n", chFile );
00387
00388 fdes = open_data( chFile, "w", AS_LOCAL_ONLY );
00389 lgErr = false;
00390
00391 (void)time(&timer);
00392 lgErr = lgErr || ( fprintf(fdes,"# this file was created by Cloudy %s (%s) on %s",
00393 t_version::Inst().chVersion,t_version::Inst().chDate,ctime(&timer)) < 0 );
00394 lgErr = lgErr || ( fprintf(fdes,"# ===========================================\n#\n") < 0 );
00395 lgErr = lgErr || ( fprintf(fdes,"%12ld # magic number opacity file\n",MAGIC_OPC) < 0 );
00396 lgErr = lgErr || ( fprintf(fdes,"%12ld # magic number rfi/mix file\n",gd.magic) < 0 );
00397 lgErr = lgErr || ( fprintf(fdes,"%12ld # magic number szd file\n",sd.magic) < 0 );
00398
00399
00400
00401 strncpy(hlp1,chFile,(size_t)(LABELSUB1+1));
00402 hlp1[LABELSUB1+1] = '\0';
00403 str = strstr_s(hlp1,"-");
00404
00405 if( str != NULL )
00406 *str = '\0';
00407
00408 strncpy(hlp2,chFile2,(size_t)(LABELSUB2+1));
00409 hlp2[LABELSUB2+1] = '\0';
00410 str = strstr_s(hlp2,"-");
00411
00412 if( str != NULL )
00413 *str = '\0';
00414
00415 strcpy(chGrainLabel," ");
00416 if( sd.nPart > 1 )
00417 {
00418 hlp1[LABELSUB1] = '\0';
00419 hlp2[LABELSUB2] = '\0';
00420 strcat(strcat(strcat(strcat(chGrainLabel,hlp1),"-"),hlp2),"xx");
00421 lgErr = lgErr || ( fprintf(fdes,"%-12.12s # grain type label, xx will be replaced by bin no.\n",
00422 chGrainLabel) < 0 );
00423 }
00424 else
00425 {
00426 strcat(strcat(strcat(chGrainLabel,hlp1),"-"),hlp2);
00427 lgErr = lgErr || ( fprintf(fdes,"%-12.12s # grain type label\n", chGrainLabel) < 0 );
00428 }
00429
00430 lgErr = lgErr || ( fprintf(fdes,"%.6e # specific weight (g/cm^3)\n",gd.rho) < 0 );
00431 lgErr = lgErr || ( fprintf(fdes,"%.6e # molecular weight of grain molecule (amu)\n",gd.mol_weight) < 0 );
00432 lgErr = lgErr || ( fprintf(fdes,"%.6e # average molecular weight per atom (amu)\n", gd.atom_weight) < 0 );
00433 lgErr = lgErr || ( fprintf(fdes,"%.6e # abundance of grain molecule at max depletion\n",gd.abun) < 0 );
00434 lgErr = lgErr || ( fprintf(fdes,"%.6e # depletion efficiency\n",gd.depl) < 0 );
00435 lgErr = lgErr || ( fprintf(fdes,"%.6e # average grain radius <a^3>/<a^2>, full size distr (cm)\n",
00436 3.*sd.vol/sd.area) < 0 );
00437 lgErr = lgErr || ( fprintf(fdes,"%.6e # average grain surface area <4pi*a^2>, full size distr (cm^2)\n",
00438 sd.area) < 0 );
00439 lgErr = lgErr || ( fprintf(fdes,"%.6e # average grain volume <4/3pi*a^3>, full size distr (cm^3)\n",
00440 sd.vol) < 0 );
00441 lgErr = lgErr || ( fprintf(fdes,"%.6e # total grain radius Int(a) per H, full size distr (cm/H)\n",
00442 sd.radius/gd.norm) < 0 );
00443 lgErr = lgErr || ( fprintf(fdes,"%.6e # total grain area Int(4pi*a^2) per H, full size distr (cm^2/H)\n",
00444 sd.area/gd.norm) < 0 );
00445 lgErr = lgErr || ( fprintf(fdes,"%.6e # total grain vol Int(4/3pi*a^3) per H, full size distr (cm^3/H)\n",
00446 sd.vol/gd.norm) < 0 );
00447 lgErr = lgErr || ( fprintf(fdes,"%.6e # work function (Ryd)\n",gd.work) < 0 );
00448 lgErr = lgErr || ( fprintf(fdes,"%.6e # gap between valence and conduction band (Ryd)\n",gd.bandgap) < 0 );
00449 lgErr = lgErr || ( fprintf(fdes,"%.6e # efficiency of thermionic emission\n",gd.therm_eff) < 0 );
00450 lgErr = lgErr || ( fprintf(fdes,"%.6e # sublimation temperature (K)\n",gd.subl_temp) < 0 );
00451 lgErr = lgErr || ( fprintf(fdes,"%12d # material type, 1=carbonaceous, 2=silicate\n",gd.matType) < 0 );
00452 lgErr = lgErr || ( fprintf(fdes,"#\n# abundances of constituent elements rel. to hydrogen\n#\n") < 0 );
00453
00454 for( nelem=0; nelem < LIMELM; nelem++ )
00455 {
00456 lgErr = lgErr || ( fprintf(fdes,"%.6e # %s\n",gd.elmAbun[nelem],
00457 elementnames.chElementSym[nelem]) < 0 );
00458 }
00459
00460 if( sd.sdCase != SD_SINGLE_SIZE && sd.sdCase != SD_NR_CARBON )
00461 {
00462 lgErr = lgErr || ( fprintf(fdes,"#\n# entire size distribution, amin=%.5f amax=%.5f micron\n",
00463 sd.lim[ipBLo],sd.lim[ipBHi]) < 0 );
00464 lgErr = lgErr || ( fprintf(fdes,"#\n%.6e # ratio a_max/a_min in each size bin\n",
00465 pow(sd.lim[ipBHi]/sd.lim[ipBLo],1./(double)sd.nPart) ) < 0 );
00466 lgErr = lgErr || ( fprintf(fdes,"#\n# size distribution function\n#\n") < 0 );
00467 lgErr = lgErr || ( fprintf(fdes,"%12ld # number of table entries\n#\n",NPTS_TABLE+1) < 0 );
00468 lgErr = lgErr || ( fprintf(fdes,"# ============================\n") < 0 );
00469 lgErr = lgErr || ( fprintf(fdes,"# size (micr) a^4*dN/da (cm^3/H)\n#\n") < 0 );
00470 for( i=0; i <= NPTS_TABLE; i++ )
00471 {
00472 double radius, a4dNda;
00473 radius = sd.lim[ipBLo]*exp((double)i/(double)NPTS_TABLE*log(sd.lim[ipBHi]/sd.lim[ipBLo]));
00474 radius = max(min(radius,sd.lim[ipBHi]),sd.lim[ipBLo]);
00475 a4dNda = POW4(radius)*size_distr(radius,&sd)/gd.norm*1.e-12/sd.unity;
00476 lgErr = lgErr || ( fprintf(fdes,"%.6e %.6e\n",radius,a4dNda) < 0 );
00477 }
00478 }
00479 else
00480 {
00481 lgErr = lgErr || ( fprintf(fdes,"#\n") < 0 );
00482 lgErr = lgErr || ( fprintf(fdes,"%.6e # a_max/a_min = 1 for single sized grain\n", 1. ) < 0 );
00483 lgErr = lgErr || ( fprintf(fdes,"%12ld # no size distribution table\n",0L) < 0 );
00484 }
00485
00486 union {
00487 double x;
00488 uint32 i[2];
00489 } u;
00490
00491 lgErr = lgErr || ( fprintf(fdes,"#\n") < 0 );
00492 lgErr = lgErr || ( fprintf(fdes,"%s # check 1\n",continuum.mesh_md5sum.c_str()) < 0 );
00493 u.x = double(rfield.emm);
00494 if( cpu.big_endian() )
00495 lgErr = lgErr || ( fprintf(fdes,"%23.8x %8.8x # check 2\n",u.i[0],u.i[1]) < 0 );
00496 else
00497 lgErr = lgErr || ( fprintf(fdes,"%23.8x %8.8x # check 2\n",u.i[1],u.i[0]) < 0 );
00498 u.x = double(rfield.egamry);
00499 if( cpu.big_endian() )
00500 lgErr = lgErr || ( fprintf(fdes,"%23.8x %8.8x # check 3\n",u.i[0],u.i[1]) < 0 );
00501 else
00502 lgErr = lgErr || ( fprintf(fdes,"%23.8x %8.8x # check 3\n",u.i[1],u.i[0]) < 0 );
00503 u.x = continuum.ResolutionScaleFactor;
00504 if( cpu.big_endian() )
00505 lgErr = lgErr || ( fprintf(fdes,"%23.8x %8.8x # check 4\n",u.i[0],u.i[1]) < 0 );
00506 else
00507 lgErr = lgErr || ( fprintf(fdes,"%23.8x %8.8x # check 4\n",u.i[1],u.i[0]) < 0 );
00508 lgErr = lgErr || ( fprintf(fdes,"%32ld # rfield.nupper\n",rfield.nupper) < 0 );
00509 lgErr = lgErr || ( fprintf(fdes,"%32ld # number of size distr. bins\n#\n",sd.nPart) < 0 );
00510
00511 if( gd.rfiType == OPC_PAH1 )
00512 {
00513 gd2.clear();
00514 mie_read_rfi("graphite.rfi",&gd2);
00515 }
00516 else if( gd.rfiType == OPC_PAH2N || gd.rfiType == OPC_PAH2C )
00517 {
00518 gd2.clear();
00519 mie_read_rfi("gdraine.rfi",&gd2);
00520 }
00521
00522 vector<int> ErrorIndex( rfield.nupper );
00523
00524 for( p=0; p < sd.nPart; p++ )
00525 {
00526 sd.cPart = p;
00527
00528 mie_auxiliary(&sd,&gd,"step");
00529
00530 if( sd.nPart > 1 )
00531 {
00532
00533
00534
00535
00536
00537 double frac = sd.unity_bin/sd.unity;
00538 volfrac = sd.vol*frac/volnorm;
00539 fprintf( ioQQQ, " Starting size bin %ld, amin=%.5f amax=%.5f micron\n",
00540 p+1,sd.clim[ipBLo],sd.clim[ipBHi] );
00541 lgErr = lgErr || ( fprintf(fdes,"# size bin %ld, amin=%.5f amax=%.5f micron\n",
00542 p+1,sd.clim[ipBLo],sd.clim[ipBHi]) < 0 );
00543 lgErr = lgErr || ( fprintf(fdes,"%.6e # average grain ",3.*sd.vol/sd.area) < 0 );
00544 lgErr = lgErr || ( fprintf(fdes,"radius <a^3>/<a^2>, this bin (cm)\n") < 0 );
00545 lgErr = lgErr || ( fprintf(fdes,"%.6e # average ",sd.area) < 0 );
00546 lgErr = lgErr || ( fprintf(fdes,"grain area <4pi*a^2>, this bin (cm^2)\n") < 0 );
00547 lgErr = lgErr || ( fprintf(fdes,"%.6e # average ",sd.vol) < 0 );
00548 lgErr = lgErr || ( fprintf(fdes,"grain volume <4/3pi*a^3>, this bin (cm^3)\n") < 0 );
00549 lgErr = lgErr || ( fprintf(fdes,"%.6e # total grain ",sd.radius*frac/gd.norm) < 0 );
00550 lgErr = lgErr || ( fprintf(fdes,"radius Int(a) per H, this bin (cm/H)\n") < 0 );
00551 lgErr = lgErr || ( fprintf(fdes,"%.6e # total grain area ",sd.area*frac/gd.norm) < 0 );
00552 lgErr = lgErr || ( fprintf(fdes,"Int(4pi*a^2) per H, this bin (cm^2/H)\n") < 0 );
00553 lgErr = lgErr || ( fprintf(fdes,"%.6e # total grain volume ",sd.vol*frac/gd.norm) < 0 );
00554 lgErr = lgErr || ( fprintf(fdes,"Int(4/3pi*a^3) per H, this bin (cm^3/H)\n#\n") < 0 );
00555 }
00556
00557 lgErrorOccurred = false;
00558
00559
00560 for( i=0; i < rfield.nupper; i++ )
00561 {
00562 wavlen = WAVNRYD/rfield.anu[i]*1.e4;
00563
00564 ErrorIndex[i] = 0;
00565 acs_abs[p][i] = 0.;
00566 acs_sct[p][i] = 0.;
00567 a1g[p][i] = 0.;
00568
00569 switch( gd.rfiType )
00570 {
00571 case RFI_TABLE:
00572 for( gd.cAxis=0; gd.cAxis < gd.nAxes; gd.cAxis++ )
00573 {
00574 mie_cs_size_distr(wavlen,&sd,&gd,mie_cs,&cs_abs,&cs_sct,&cosb,&Error);
00575 ErrorIndex[i] = max(ErrorIndex[i],Error);
00576 acs_abs[p][i] += cs_abs*gd.wt[gd.cAxis];
00577 acs_sct[p][i] += cs_sct*gd.wt[gd.cAxis];
00578 a1g[p][i] += cs_sct*(1.-cosb)*gd.wt[gd.cAxis];
00579 }
00580 lgErrorOccurred = mie_auxiliary2(ErrorIndex,acs_abs,acs_sct,a1g,p,i);
00581 break;
00582 case OPC_TABLE:
00583 gd.cAxis = 0;
00584 mie_cs_size_distr(wavlen,&sd,&gd,tbl_fun,&cs_abs,&cs_sct,&cosb,&Error);
00585 ErrorIndex[i] = min(Error,2);
00586 lgErrorOccurred = lgErrorOccurred || ( Error > 0 );
00587 acs_abs[p][i] = cs_abs*gd.norm;
00588 acs_sct[p][i] = cs_sct*gd.norm;
00589 a1g[p][i] = 1.-cosb;
00590 break;
00591 case OPC_GREY:
00592 ErrorIndex[i] = 0;
00593 acs_abs[p][i] = 1.3121e-23*volfrac*gd.norm;
00594 acs_sct[p][i] = 2.6242e-23*volfrac*gd.norm;
00595 a1g[p][i] = 1.;
00596 break;
00597 case OPC_PAH1:
00598 gd.cAxis = 0;
00599 for( gd2.cAxis=0; gd2.cAxis < gd2.nAxes; gd2.cAxis++ )
00600 {
00601 mie_cs_size_distr(wavlen,&sd,&gd,car1_fun,&cs_abs,&cs_sct,&cosb,&Error);
00602 ErrorIndex[i] = max(ErrorIndex[i],Error);
00603 acs_abs[p][i] += cs_abs*gd2.wt[gd2.cAxis];
00604 acs_sct[p][i] += 0.1*cs_abs*gd2.wt[gd2.cAxis];
00605 a1g[p][i] += 0.1*cs_abs*1.*gd2.wt[gd2.cAxis];
00606 }
00607 lgErrorOccurred = mie_auxiliary2(ErrorIndex,acs_abs,acs_sct,a1g,p,i);
00608 break;
00609 case OPC_PAH2N:
00610 case OPC_PAH2C:
00611 gd.cAxis = 0;
00612
00613 gd.charge = ( gd.rfiType == OPC_PAH2N ) ? 0 : 1;
00614 for( gd2.cAxis=0; gd2.cAxis < gd2.nAxes; gd2.cAxis++ )
00615 {
00616 mie_cs_size_distr(wavlen,&sd,&gd,ld01_fun,&cs_abs,&cs_sct,&cosb,&Error);
00617 ErrorIndex[i] = max(ErrorIndex[i],Error);
00618 acs_abs[p][i] += cs_abs*gd2.wt[gd2.cAxis];
00619 acs_sct[p][i] += 0.1*cs_abs*gd2.wt[gd2.cAxis];
00620 a1g[p][i] += 0.1*cs_abs*1.*gd2.wt[gd2.cAxis];
00621 }
00622 lgErrorOccurred = mie_auxiliary2(ErrorIndex,acs_abs,acs_sct,a1g,p,i);
00623 break;
00624 default:
00625 fprintf( ioQQQ, " Insanity in mie_write_opc\n" );
00626 ShowMe();
00627 cdEXIT(EXIT_FAILURE);
00628 }
00629 }
00630
00631
00632 if( lgErrorOccurred )
00633 {
00634 strcpy(chString,"absorption cs");
00635 mie_repair(chString,rfield.nupper,2,0,rfield.anu,&acs_abs[p][0],ErrorIndex,false,&lgWarning);
00636 strcpy(chString,"scattering cs");
00637 mie_repair(chString,rfield.nupper,2,1,rfield.anu,&acs_sct[p][0],ErrorIndex,false,&lgWarning);
00638 strcpy(chString,"asymmetry parameter");
00639 mie_repair(chString,rfield.nupper,1,1,rfield.anu,&a1g[p][0],ErrorIndex,true,&lgWarning);
00640 }
00641
00642 for( i=0; i < rfield.nupper; i++ )
00643 {
00644 acs_abs[p][i] /= gd.norm;
00645
00646
00647 acs_sct[p][i] /= gd.norm;
00648 }
00649 }
00650
00651 lgErr = lgErr || ( fprintf(fdes,"#\n") < 0 );
00652 lgErr = lgErr || ( fprintf(fdes,"# ===========================================\n") < 0 );
00653 lgErr = lgErr || ( fprintf(fdes,"# anu (Ryd) abs_cs_01 (cm^2/H) abs_cs_02.....\n#\n") < 0 );
00654
00655 for( i=0; i < rfield.nupper; i++ )
00656 {
00657 lgErr = lgErr || ( fprintf(fdes,"%.6e ",rfield.anu[i]) < 0 );
00658 for( p=0; p < sd.nPart; p++ )
00659 {
00660 lgErr = lgErr || ( fprintf(fdes,"%.6e ",acs_abs[p][i]) < 0 );
00661 }
00662 lgErr = lgErr || ( fprintf(fdes,"\n") < 0 );
00663 }
00664
00665 lgErr = lgErr || ( fprintf(fdes,"#\n") < 0 );
00666 lgErr = lgErr || ( fprintf(fdes,"# ===========================================\n") < 0 );
00667 lgErr = lgErr || ( fprintf(fdes,"# anu (Ryd) sct_cs_01 (cm^2/H) sct_cs_02.....\n#\n") < 0 );
00668
00669 for( i=0; i < rfield.nupper; i++ )
00670 {
00671 lgErr = lgErr || ( fprintf(fdes,"%.6e ",rfield.anu[i]) < 0 );
00672 for( p=0; p < sd.nPart; p++ )
00673 {
00674 lgErr = lgErr || ( fprintf(fdes,"%.6e ",acs_sct[p][i]) < 0 );
00675 }
00676 lgErr = lgErr || ( fprintf(fdes,"\n") < 0 );
00677 }
00678
00679 lgErr = lgErr || ( fprintf(fdes,"#\n") < 0 );
00680 lgErr = lgErr || ( fprintf(fdes,"# ===========================================\n") < 0 );
00681 lgErr = lgErr || ( fprintf(fdes,"# anu (Ryd) (1-g)_bin_01 (1-g)_bin_02.....\n#\n") < 0 );
00682
00683 for( i=0; i < rfield.nupper; i++ )
00684 {
00685 lgErr = lgErr || ( fprintf(fdes,"%.6e ",rfield.anu[i]) < 0 );
00686 for( p=0; p < sd.nPart; p++ )
00687 {
00688 lgErr = lgErr || ( fprintf(fdes,"%.6e ",a1g[p][i]) < 0 );
00689 }
00690 lgErr = lgErr || ( fprintf(fdes,"\n") < 0 );
00691 }
00692
00693 fprintf( ioQQQ, " Starting calculation of inverse attenuation length\n" );
00694 strcpy(chString,"inverse attenuation length");
00695 if( gd.rfiType != RFI_TABLE )
00696 {
00697
00698 gd2.clear();
00699 ial_type icase = gv.which_ial[gd.matType];
00700 switch( icase )
00701 {
00702 case IAL_CAR:
00703 mie_read_rfi("graphite.rfi",&gd2);
00704 mie_calc_ial(&gd2,rfield.nupper,inv_att_len,chString,&lgWarning);
00705 break;
00706 case IAL_SIL:
00707 mie_read_rfi("silicate.rfi",&gd2);
00708 mie_calc_ial(&gd2,rfield.nupper,inv_att_len,chString,&lgWarning);
00709 break;
00710 default:
00711 fprintf( ioQQQ, " mie_write_opc detected unknown ial type: %d\n" , icase );
00712 cdEXIT(EXIT_FAILURE);
00713 }
00714 }
00715 else
00716 {
00717 mie_calc_ial(&gd,rfield.nupper,inv_att_len,chString,&lgWarning);
00718 }
00719
00720 lgErr = lgErr || ( fprintf(fdes,"#\n") < 0 );
00721 lgErr = lgErr || ( fprintf(fdes,"# ===========================================\n") < 0 );
00722 lgErr = lgErr || ( fprintf(fdes,"# anu (Ryd) inverse attenuation length (cm^-1)\n#\n") < 0 );
00723
00724 for( i=0; i < rfield.nupper; i++ )
00725 {
00726 lgErr = lgErr || ( fprintf(fdes,"%.6e %.6e\n",rfield.anu[i],inv_att_len[i]) < 0 );
00727 }
00728
00729 fclose(fdes);
00730
00731 if( lgErr )
00732 {
00733 fprintf( ioQQQ, "\n Error writing file: %s\n", chFile );
00734 if( remove(chFile) == 0 )
00735 {
00736 fprintf( ioQQQ, " The file has been removed\n" );
00737 cdEXIT(EXIT_FAILURE);
00738 }
00739 }
00740 else
00741 {
00742 fprintf( ioQQQ, "\n Opacity file %s written succesfully\n\n", chFile );
00743 if( lgWarning )
00744 {
00745 fprintf( ioQQQ, "\n !!! Warnings were detected !!!\n\n" );
00746 }
00747 }
00748 return;
00749 }
00750
00751
00752 STATIC void mie_auxiliary( sd_data *sd,
00753 const grain_data *gd,
00754 const char *auxCase)
00755 {
00756 double amin,
00757 amax,
00758 delta,
00759 oldvol,
00760 step;
00761
00762
00763 const double TOLER = 1.e-3;
00764
00765 DEBUG_ENTRY( "mie_auxiliary()" );
00766 if( strcmp(auxCase,"init") == 0 )
00767 {
00768 double mass, radius, CpMolecule;
00769
00770
00771
00772 sd->nmul = 1;
00773
00774
00775 switch( sd->sdCase )
00776 {
00777 case SD_SINGLE_SIZE:
00778 sd->radius = sd->a[ipSize]*1.e-4;
00779 sd->area = 4.*PI*pow2(sd->a[ipSize])*1.e-8;
00780 sd->vol = 4./3.*PI*pow3(sd->a[ipSize])*1.e-12;
00781 break;
00782 case SD_NR_CARBON:
00783 if( gd->elmAbun[ipCARBON] == 0. )
00784 {
00785 fprintf( ioQQQ, "\n This size distribution can only be combined with"
00786 " carbonaceous material, bailing out...\n" );
00787 cdEXIT(EXIT_FAILURE);
00788 }
00789
00790 CpMolecule = gd->elmAbun[ipCARBON]/(gd->abun*gd->depl);
00791
00792 mass = (double)sd->nCarbon/CpMolecule*gd->mol_weight*ATOMIC_MASS_UNIT;
00793 radius = pow(3.*mass/(PI4*gd->rho),1./3.);
00794 sd->a[ipSize] = radius*1.e4;
00795 sd->radius = radius;
00796 sd->area = 4.*PI*pow2(radius);
00797 sd->vol = 4./3.*PI*pow3(radius);
00798 break;
00799 case SD_POWERLAW:
00800 case SD_EXP_CUTOFF1:
00801 case SD_EXP_CUTOFF2:
00802 case SD_EXP_CUTOFF3:
00803 case SD_LOG_NORMAL:
00804 case SD_LIN_NORMAL:
00805 case SD_TABLE:
00806
00807
00808 amin = sd->lgLogScale ? log(sd->lim[ipBLo]) : sd->lim[ipBLo];
00809 amax = sd->lgLogScale ? log(sd->lim[ipBHi]) : sd->lim[ipBHi];
00810
00811 sd->clim[ipBLo] = sd->lim[ipBLo];
00812 sd->clim[ipBHi] = sd->lim[ipBHi];
00813
00814
00815 oldvol= 0.;
00816 do
00817 {
00818 sd->nmul *= 2;
00819 mie_integrate(sd,amin,amax,&sd->unity);
00820 delta = fabs(sd->vol-oldvol)/sd->vol;
00821 oldvol = sd->vol;
00822 } while( sd->nmul <= 1024 && delta > TOLER );
00823
00824 if( delta > TOLER )
00825 {
00826 fprintf( ioQQQ, " could not converge integration of size distribution\n" );
00827 cdEXIT(EXIT_FAILURE);
00828 }
00829
00830
00831
00832 sd->nmul /= 2;
00833 mie_integrate(sd,amin,amax,&sd->unity);
00834 break;
00835 case SD_ILLEGAL:
00836 default:
00837 fprintf( ioQQQ, " insane case for grain size distribution: %d\n" , sd->sdCase );
00838 ShowMe();
00839 cdEXIT(EXIT_FAILURE);
00840 }
00841 }
00842 else if( strcmp(auxCase,"step") == 0 )
00843 {
00844
00845 switch( sd->sdCase )
00846 {
00847 case SD_SINGLE_SIZE:
00848 case SD_NR_CARBON:
00849 break;
00850 case SD_POWERLAW:
00851 case SD_EXP_CUTOFF1:
00852 case SD_EXP_CUTOFF2:
00853 case SD_EXP_CUTOFF3:
00854 case SD_LOG_NORMAL:
00855 case SD_LIN_NORMAL:
00856 case SD_TABLE:
00857 amin = sd->lgLogScale ? log(sd->lim[ipBLo]) : sd->lim[ipBLo];
00858 amax = sd->lgLogScale ? log(sd->lim[ipBHi]) : sd->lim[ipBHi];
00859 step = (amax - amin)/(double)sd->nPart;
00860 amin = amin + (double)sd->cPart*step;
00861 amax = min(amax,amin + step);
00862
00863 sd->clim[ipBLo] = sd->lgLogScale ? exp(amin) : amin;
00864 sd->clim[ipBHi] = sd->lgLogScale ? exp(amax) : amax;
00865
00866 mie_integrate(sd,amin,amax,&sd->unity_bin);
00867
00868 break;
00869 case SD_ILLEGAL:
00870 default:
00871 fprintf( ioQQQ, " insane case for grain size distribution: %d\n" , sd->sdCase );
00872 ShowMe();
00873 cdEXIT(EXIT_FAILURE);
00874 }
00875 }
00876 else
00877 {
00878 fprintf( ioQQQ, " mie_auxiliary called with insane argument: %s\n", auxCase );
00879 ShowMe();
00880 cdEXIT(EXIT_FAILURE);
00881 }
00882 return;
00883 }
00884
00885
00886 STATIC bool mie_auxiliary2(vector<int>& ErrorIndex,
00887 multi_arr<double,2>& acs_abs,
00888 multi_arr<double,2>& acs_sct,
00889 multi_arr<double,2>& a1g,
00890 long p,
00891 long i)
00892 {
00893 DEBUG_ENTRY( "mie_auxiliary2()" );
00894
00895 bool lgErrorOccurred = false;
00896 if( ErrorIndex[i] > 0 )
00897 {
00898 ErrorIndex[i] = min(ErrorIndex[i],2);
00899 lgErrorOccurred = true;
00900 }
00901
00902 switch( ErrorIndex[i] )
00903 {
00904
00905 case 2:
00906 acs_abs[p][i] = 0.;
00907 acs_sct[p][i] = 0.;
00908
00909
00910 case 1:
00911 a1g[p][i] = 0.;
00912 break;
00913
00914 case 0:
00915 a1g[p][i] /= acs_sct[p][i];
00916 break;
00917 default:
00918 fprintf( ioQQQ, " Insane value for ErrorIndex: %d\n", ErrorIndex[i] );
00919 ShowMe();
00920 cdEXIT(EXIT_FAILURE);
00921 }
00922
00923
00924 if( ErrorIndex[i] < 2 )
00925 ASSERT( acs_abs[p][i] > 0. && acs_sct[p][i] > 0. );
00926 if( ErrorIndex[i] < 1 )
00927 ASSERT( a1g[p][i] > 0. );
00928
00929 return lgErrorOccurred;
00930 }
00931
00932
00933 STATIC void mie_integrate( sd_data *sd,
00934 double amin,
00935 double amax,
00936 double *normalization)
00937 {
00938 long int j;
00939 double unity;
00940
00941 DEBUG_ENTRY( "mie_integrate()" );
00942
00943
00944
00945 sd->nn = sd->nmul*((long)(2.*log(sd->clim[ipBHi]/sd->clim[ipBLo])) + 1);
00946 sd->nn = min(max(sd->nn,2*sd->nmul),4096);
00947 sd->xx.resize(sd->nn);
00948 sd->aa.resize(sd->nn);
00949 sd->rr.resize(sd->nn);
00950 sd->ww.resize(sd->nn);
00951 gauss_legendre(sd->nn,sd->xx,sd->aa);
00952 gauss_init(sd->nn,amin,amax,sd->xx,sd->aa,sd->rr,sd->ww);
00953
00954
00955 unity = 0.;
00956 sd->radius = 0.;
00957 sd->area = 0.;
00958 sd->vol = 0.;
00959
00960 for( j=0; j < sd->nn; j++ )
00961 {
00962 double weight;
00963
00964
00965 if( sd->lgLogScale )
00966 {
00967 sd->rr[j] = exp(sd->rr[j]);
00968 sd->ww[j] *= sd->rr[j];
00969 }
00970 weight = sd->ww[j]*size_distr(sd->rr[j],sd);
00971 unity += weight;
00972 sd->radius += weight*sd->rr[j];
00973 sd->area += weight*pow2(sd->rr[j]);
00974 sd->vol += weight*pow3(sd->rr[j]);
00975 }
00976 *normalization = unity;
00977 sd->radius *= 1.e-4/unity;
00978 sd->area *= 4.*PI*1.e-8/unity;
00979 sd->vol *= 4./3.*PI*1.e-12/unity;
00980 return;
00981 }
00982
00983
00984 void mie_read_opc(const char *chFile,
00985 const GrainPar& gp)
00986 {
00987 int res,
00988 lgDefaultQHeat;
00989 long int dl,
00990 help,
00991 i,
00992 nelem,
00993 j,
00994 magic,
00995 nbin,
00996 nup;
00997 size_t nd,
00998 nd2;
00999 realnum RefAbund[LIMELM],
01000 VolTotal;
01001 double anu;
01002 double RadiusRatio;
01003 char chLine[FILENAME_PATH_LENGTH_2],
01004 md5sum[NMD5+1],
01005 *str;
01006 FILE *io2;
01007
01008
01009
01010 const double RATIO_MAX = pow(100.,1./3.);
01011
01012 DEBUG_ENTRY( "mie_read_opc()" );
01013
01014 io2 = open_data( chFile, "r", AS_DATA_LOCAL );
01015
01016
01017 strcpy( chLine, " " );
01018 if( strlen(chFile) <= 40 )
01019 {
01020 strncpy( &chLine[0], chFile, strlen(chFile) );
01021 }
01022 else
01023 {
01024 strncpy( &chLine[0], chFile, 37 );
01025 strncpy( &chLine[37], "...", 3 );
01026 }
01027 if( called.lgTalk )
01028 fprintf( ioQQQ, " * >>>> mie_read_opc reading file -- %40s <<<< *\n", chLine );
01029
01030
01031 for( size_t p=0; p < gv.ReadRecord.size(); ++p )
01032 {
01033 if( gv.ReadRecord[p] == chFile )
01034 {
01035 fprintf( ioQQQ, " File %s has already been read before, was this intended ?\n", chFile );
01036 break;
01037 }
01038 }
01039 gv.ReadRecord.push_back( chFile );
01040
01041
01042 gv.bin.push_back( new GrainBin );
01043 nd = gv.bin.size()-1;
01044
01045 dl = 0;
01046
01047
01048 mie_next_data(chFile,io2,chLine,&dl);
01049 mie_read_long(chFile,chLine,&magic,true,dl);
01050 if( magic != MAGIC_OPC )
01051 {
01052 fprintf( ioQQQ, " Opacity file %s has obsolete magic number\n",chFile );
01053 fprintf( ioQQQ, " I found magic number %ld, but expected %ld on line #%ld\n",
01054 magic,MAGIC_OPC,dl );
01055 fprintf( ioQQQ, " Please recompile this file\n" );
01056 cdEXIT(EXIT_FAILURE);
01057 }
01058
01059
01060 mie_next_data(chFile,io2,chLine,&dl);
01061 mie_read_long(chFile,chLine,&magic,true,dl);
01062
01063 mie_next_data(chFile,io2,chLine,&dl);
01064 mie_read_long(chFile,chLine,&magic,true,dl);
01065
01066
01067
01068
01069
01070
01071 gv.bin[nd]->dstfactor = (realnum)gp.dep;
01072
01073
01074 mie_next_data(chFile,io2,chLine,&dl);
01075 strncpy(gv.bin[nd]->chDstLab,chLine,(size_t)LABELSIZE);
01076 gv.bin[nd]->chDstLab[LABELSIZE] = '\0';
01077
01078
01079 mie_next_data(chFile,io2,chLine,&dl);
01080 mie_read_realnum(chFile,chLine,&gv.bin[nd]->dustp[0],true,dl);
01081
01082 gv.bin[nd]->eec = pow((double)gv.bin[nd]->dustp[0],-0.85);
01083
01084
01085 mie_next_data(chFile,io2,chLine,&dl);
01086 mie_read_realnum(chFile,chLine,&gv.bin[nd]->dustp[1],true,dl);
01087
01088
01089 mie_next_data(chFile,io2,chLine,&dl);
01090 mie_read_realnum(chFile,chLine,&gv.bin[nd]->atomWeight,true,dl);
01091
01092
01093 mie_next_data(chFile,io2,chLine,&dl);
01094 mie_read_realnum(chFile,chLine,&gv.bin[nd]->dustp[2],true,dl);
01095
01096
01097 mie_next_data(chFile,io2,chLine,&dl);
01098 mie_read_realnum(chFile,chLine,&gv.bin[nd]->dustp[3],true,dl);
01099
01100
01101 gv.bin[nd]->dustp[4] = 1.;
01102
01103
01104 mie_next_data(chFile,io2,chLine,&dl);
01105 mie_read_realnum(chFile,chLine,&gv.bin[nd]->AvRadius,true,dl);
01106 gv.bin[nd]->eyc = 1./gv.bin[nd]->AvRadius + 1.e7;
01107
01108
01109 mie_next_data(chFile,io2,chLine,&dl);
01110 mie_read_realnum(chFile,chLine,&gv.bin[nd]->AvArea,true,dl);
01111
01112
01113 mie_next_data(chFile,io2,chLine,&dl);
01114 mie_read_realnum(chFile,chLine,&gv.bin[nd]->AvVol,true,dl);
01115
01116
01117 mie_next_data(chFile,io2,chLine,&dl);
01118 mie_read_realnum(chFile,chLine,&gv.bin[nd]->IntRadius,true,dl);
01119
01120
01121 mie_next_data(chFile,io2,chLine,&dl);
01122 mie_read_realnum(chFile,chLine,&gv.bin[nd]->IntArea,true,dl);
01123
01124
01125 mie_next_data(chFile,io2,chLine,&dl);
01126 mie_read_realnum(chFile,chLine,&gv.bin[nd]->IntVol,true,dl);
01127
01128
01129 mie_next_data(chFile,io2,chLine,&dl);
01130 mie_read_realnum(chFile,chLine,&gv.bin[nd]->DustWorkFcn,true,dl);
01131
01132
01133 mie_next_data(chFile,io2,chLine,&dl);
01134 mie_read_realnum(chFile,chLine,&gv.bin[nd]->BandGap,false,dl);
01135
01136
01137 mie_next_data(chFile,io2,chLine,&dl);
01138 mie_read_realnum(chFile,chLine,&gv.bin[nd]->ThermEff,true,dl);
01139
01140
01141 mie_next_data(chFile,io2,chLine,&dl);
01142 mie_read_realnum(chFile,chLine,&gv.bin[nd]->Tsublimat,true,dl);
01143
01144
01145 mie_next_data(chFile,io2,chLine,&dl);
01146 mie_read_long(chFile,chLine,&help,true,dl);
01147 gv.bin[nd]->matType = (mat_type)help;
01148
01149 for( nelem=0; nelem < LIMELM; nelem++ )
01150 {
01151 mie_next_data(chFile,io2,chLine,&dl);
01152 mie_read_realnum(chFile,chLine,&RefAbund[nelem],false,dl);
01153
01154 gv.bin[nd]->elmAbund[nelem] = RefAbund[nelem];
01155
01156
01157 gv.bin[nd]->AccomCoef[nelem] = 2.*gv.bin[nd]->atomWeight*dense.AtomicWeight[nelem]/
01158 pow2(gv.bin[nd]->atomWeight+dense.AtomicWeight[nelem]);
01159 }
01160
01161
01162 mie_next_data(chFile,io2,chLine,&dl);
01163 mie_read_double(chFile,chLine,&RadiusRatio,true,dl);
01164
01165 gv.bin[nd]->nDustFunc = gp.nDustFunc;
01166 lgDefaultQHeat = ( RadiusRatio < RATIO_MAX && !gp.lgGreyGrain );
01167 gv.bin[nd]->lgQHeat = ( gp.lgForbidQHeating ) ? false : ( gp.lgRequestQHeating || lgDefaultQHeat );
01168 gv.bin[nd]->cnv_H_pGR = gv.bin[nd]->AvVol/gv.bin[nd]->IntVol;
01169 gv.bin[nd]->cnv_GR_pH = 1./gv.bin[nd]->cnv_H_pGR;
01170
01171
01172 gv.bin[nd]->Capacity = PI4*ELECTRIC_CONST*gv.bin[nd]->IntRadius/100.*gv.bin[nd]->cnv_H_pGR;
01173
01174
01175 mie_next_data(chFile,io2,chLine,&dl);
01176 mie_read_long(chFile,chLine,&nup,false,dl);
01177 for( i=0; i < nup; i++ )
01178 mie_next_data(chFile,io2,chLine,&dl);
01179
01180
01181 mie_next_data(chFile,io2,chLine,&dl);
01182 mie_read_word(chLine,md5sum,sizeof(md5sum),false);
01183
01184 union {
01185 double x;
01186 uint32 i[2];
01187 } u;
01188 double mesh_lo, mesh_hi;
01189
01190
01191 mie_next_data(chFile,io2,chLine,&dl);
01192 if( cpu.big_endian() )
01193 sscanf( chLine, "%x %x", &u.i[0], &u.i[1] );
01194 else
01195 sscanf( chLine, "%x %x", &u.i[1], &u.i[0] );
01196 mesh_lo = u.x;
01197
01198
01199 mie_next_data(chFile,io2,chLine,&dl);
01200 if( cpu.big_endian() )
01201 sscanf( chLine, "%x %x", &u.i[0], &u.i[1] );
01202 else
01203 sscanf( chLine, "%x %x", &u.i[1], &u.i[0] );
01204 mesh_hi = u.x;
01205
01206 if( strncmp( md5sum, continuum.mesh_md5sum.c_str(), NMD5 ) != 0 ||
01207 !fp_equal( mesh_lo, double(rfield.emm) ) ||
01208 !fp_equal( mesh_hi, double(rfield.egamry) ) )
01209 {
01210 fprintf( ioQQQ, " Opacity file %s has an incompatible energy grid.\n", chFile );
01211 fprintf( ioQQQ, " Please recompile this file using the COMPILE GRAINS command.\n" );
01212 cdEXIT(EXIT_FAILURE);
01213 }
01214
01215
01216 mie_next_data(chFile,io2,chLine,&dl);
01217 if( cpu.big_endian() )
01218 sscanf( chLine, "%x %x", &u.i[0], &u.i[1] );
01219 else
01220 sscanf( chLine, "%x %x", &u.i[1], &u.i[0] );
01221
01222 gv.bin[nd]->RSFCheck = u.x;
01223
01224
01225 mie_next_data(chFile,io2,chLine,&dl);
01226 mie_read_long(chFile,chLine,&nup,true,dl);
01227
01228
01229 mie_next_data(chFile,io2,chLine,&dl);
01230 mie_read_long(chFile,chLine,&nbin,true,dl);
01231
01232
01233 if( nbin > 1 )
01234 {
01235
01236 VolTotal = gv.bin[nd]->IntVol;
01237
01238 for( i=0; i < nbin; i++ )
01239 {
01240 if( i == 0 )
01241 nd2 = nd;
01242 else
01243 {
01244
01245 gv.bin.push_back( new GrainBin );
01246 nd2 = gv.bin.size()-1;
01247
01248
01249 *gv.bin[nd2] = *gv.bin[nd];
01250 }
01251
01252
01253
01254
01255 mie_next_data(chFile,io2,chLine,&dl);
01256 mie_read_realnum(chFile,chLine,&gv.bin[nd2]->AvRadius,true,dl);
01257
01258
01259 mie_next_data(chFile,io2,chLine,&dl);
01260 mie_read_realnum(chFile,chLine,&gv.bin[nd2]->AvArea,true,dl);
01261
01262
01263 mie_next_data(chFile,io2,chLine,&dl);
01264 mie_read_realnum(chFile,chLine,&gv.bin[nd2]->AvVol,true,dl);
01265
01266
01267 mie_next_data(chFile,io2,chLine,&dl);
01268 mie_read_realnum(chFile,chLine,&gv.bin[nd2]->IntRadius,true,dl);
01269
01270
01271 mie_next_data(chFile,io2,chLine,&dl);
01272 mie_read_realnum(chFile,chLine,&gv.bin[nd2]->IntArea,true,dl);
01273
01274
01275 mie_next_data(chFile,io2,chLine,&dl);
01276 mie_read_realnum(chFile,chLine,&gv.bin[nd2]->IntVol,true,dl);
01277
01278 gv.bin[nd2]->cnv_H_pGR = gv.bin[nd2]->AvVol/gv.bin[nd2]->IntVol;
01279 gv.bin[nd2]->cnv_GR_pH = 1./gv.bin[nd2]->cnv_H_pGR;
01280
01281
01282 gv.bin[nd2]->Capacity =
01283 PI4*ELECTRIC_CONST*gv.bin[nd2]->IntRadius/100.*gv.bin[nd2]->cnv_H_pGR;
01284
01285
01286
01287
01288 gv.bin[nd2]->dustp[4] = gv.bin[nd2]->IntVol/VolTotal;
01289 for( nelem=0; nelem < LIMELM; nelem++ )
01290 gv.bin[nd2]->elmAbund[nelem] = RefAbund[nelem]*gv.bin[nd2]->dustp[4];
01291 }
01292
01293
01294 for( i=0; i < nbin; i++ )
01295 {
01296 nd2 = nd + i;
01297
01298 str = strstr_s(gv.bin[nd2]->chDstLab,"xx");
01299 if( str != NULL )
01300 sprintf(str,"%02ld",i+1);
01301 }
01302 }
01303
01304
01305 for( i=0; i < nbin; i++ )
01306 {
01307 nd2 = nd + i;
01308 gv.bin[nd2]->dstab1.resize(nup);
01309 gv.bin[nd2]->pure_sc1.resize(nup);
01310 gv.bin[nd2]->asym.resize(nup);
01311 gv.bin[nd2]->inv_att_len.resize(nup);
01312 }
01313
01314
01315 for( i=0; i < 5; i++ )
01316 mie_next_line(chFile,io2,chLine,&dl);
01317
01318
01319 for( i=0; i < nup; i++ )
01320 {
01321
01322 if( (res = fscanf(io2,"%le",&anu)) != 1 )
01323 {
01324 fprintf( ioQQQ, " Read failed on %s\n",chFile );
01325 if( res == EOF )
01326 fprintf( ioQQQ, " EOF reached prematurely\n" );
01327 cdEXIT(EXIT_FAILURE);
01328 }
01329 for( j=0; j < nbin; j++ )
01330 {
01331 nd2 = nd + j;
01332 if( (res = fscanf(io2,"%le",&gv.bin[nd2]->dstab1[i])) != 1 )
01333 {
01334 fprintf( ioQQQ, " Read failed on %s\n",chFile );
01335 if( res == EOF )
01336 fprintf( ioQQQ, " EOF reached prematurely\n" );
01337 cdEXIT(EXIT_FAILURE);
01338 }
01339 ASSERT( gv.bin[nd2]->dstab1[i] > 0. );
01340 }
01341 }
01342
01343
01344 for( i=0; i < 5; i++ )
01345 mie_next_line(chFile,io2,chLine,&dl);
01346
01347
01348 for( i=0; i < nup; i++ )
01349 {
01350 if( (res = fscanf(io2,"%le",&anu)) != 1 )
01351 {
01352 fprintf( ioQQQ, " Read failed on %s\n",chFile );
01353 if( res == EOF )
01354 fprintf( ioQQQ, " EOF reached prematurely\n" );
01355 cdEXIT(EXIT_FAILURE);
01356 }
01357 for( j=0; j < nbin; j++ )
01358 {
01359 nd2 = nd + j;
01360 if( (res = fscanf(io2,"%le",&gv.bin[nd2]->pure_sc1[i])) != 1 )
01361 {
01362 fprintf( ioQQQ, " Read failed on %s\n",chFile );
01363 if( res == EOF )
01364 fprintf( ioQQQ, " EOF reached prematurely\n" );
01365 cdEXIT(EXIT_FAILURE);
01366 }
01367 ASSERT( gv.bin[nd2]->pure_sc1[i] > 0. );
01368 }
01369 }
01370
01371
01372 for( i=0; i < 5; i++ )
01373 mie_next_line(chFile,io2,chLine,&dl);
01374
01375
01376 for( i=0; i < nup; i++ )
01377 {
01378 if( (res = fscanf(io2,"%le",&anu)) != 1 )
01379 {
01380 fprintf( ioQQQ, " Read failed on %s\n",chFile );
01381 if( res == EOF )
01382 fprintf( ioQQQ, " EOF reached prematurely\n" );
01383 cdEXIT(EXIT_FAILURE);
01384 }
01385 for( j=0; j < nbin; j++ )
01386 {
01387 nd2 = nd + j;
01388 if( (res = fscanf(io2,"%le",&gv.bin[nd2]->asym[i])) != 1 )
01389 {
01390 fprintf( ioQQQ, " Read failed on %s\n",chFile );
01391 if( res == EOF )
01392 fprintf( ioQQQ, " EOF reached prematurely\n" );
01393 cdEXIT(EXIT_FAILURE);
01394 }
01395 ASSERT( gv.bin[nd2]->asym[i] > 0. );
01396 }
01397 }
01398
01399
01400 for( i=0; i < 5; i++ )
01401 mie_next_line(chFile,io2,chLine,&dl);
01402
01403
01404 for( i=0; i < nup; i++ )
01405 {
01406 double help;
01407 if( (res = fscanf(io2,"%le %le",&anu,&help)) != 2 )
01408 {
01409 fprintf( ioQQQ, " Read failed on %s\n",chFile );
01410 if( res == EOF )
01411 fprintf( ioQQQ, " EOF reached prematurely\n" );
01412 cdEXIT(EXIT_FAILURE);
01413 }
01414 gv.bin[nd]->inv_att_len[i] = (realnum)help;
01415 ASSERT( gv.bin[nd]->inv_att_len[i] > 0. );
01416
01417 for( j=1; j < nbin; j++ )
01418 {
01419 nd2 = nd + j;
01420 gv.bin[nd2]->inv_att_len[i] = gv.bin[nd]->inv_att_len[i];
01421 }
01422 }
01423
01424 fclose(io2);
01425 return;
01426 }
01427
01428
01429
01430
01431 STATIC void mie_cs_size_distr(double wavlen,
01432 sd_data *sd,
01433 const grain_data *gd,
01434 void(*cs_fun)(double,const sd_data*,
01435 const grain_data*,
01436 double*,double*,
01437 double*,int*),
01438 double *cs_abs,
01439 double *cs_sct,
01440 double *cosb,
01441 int *error)
01442 {
01443 bool lgADLused;
01444 long int i;
01445 double absval,
01446 g,
01447 sctval,
01448 weight;
01449
01450 DEBUG_ENTRY( "mie_cs_size_distr()" );
01451
01452
01453 ASSERT( wavlen > 0. );
01454 ASSERT( gd->cAxis >= 0 && gd->cAxis < gd->nAxes && gd->cAxis < NAX );
01455
01456 switch( sd->sdCase )
01457 {
01458 case SD_SINGLE_SIZE:
01459 case SD_NR_CARBON:
01460
01461 ASSERT( sd->a[ipSize] > 0. );
01462 sd->cSize = sd->a[ipSize];
01463 (*cs_fun)(wavlen,sd,gd,cs_abs,cs_sct,cosb,error);
01464 break;
01465 case SD_POWERLAW:
01466
01467 case SD_EXP_CUTOFF1:
01468 case SD_EXP_CUTOFF2:
01469 case SD_EXP_CUTOFF3:
01470
01471 case SD_LOG_NORMAL:
01472
01473 case SD_LIN_NORMAL:
01474
01475 case SD_TABLE:
01476
01477 ASSERT( sd->lim[ipBLo] > 0. && sd->lim[ipBHi] > 0. && sd->lim[ipBHi] > sd->lim[ipBLo] );
01478 lgADLused = false;
01479 *cs_abs = 0.;
01480 *cs_sct = 0.;
01481 *cosb = 0.;
01482 for( i=0; i < sd->nn; i++ )
01483 {
01484 sd->cSize = sd->rr[i];
01485 (*cs_fun)(wavlen,sd,gd,&absval,&sctval,&g,error);
01486 if( *error >= 2 )
01487 {
01488
01489 *cs_abs = -1.;
01490 *cs_sct = -1.;
01491 *cosb = -2.;
01492 return;
01493 }
01494 else if( *error == 1 )
01495 {
01496
01497 lgADLused = true;
01498 }
01499 weight = sd->ww[i]*size_distr(sd->rr[i],sd);
01500 *cs_abs += weight*absval;
01501 *cs_sct += weight*sctval;
01502 *cosb += weight*sctval*g;
01503 }
01504 if( lgADLused )
01505 {
01506 *error = 1;
01507 *cosb = -2.;
01508 }
01509 else
01510 {
01511 *error = 0;
01512 *cosb /= *cs_sct;
01513 }
01514 *cs_abs /= sd->unity;
01515 *cs_sct /= sd->unity;
01516 break;
01517 case SD_ILLEGAL:
01518 default:
01519 fprintf( ioQQQ, " insane case for grain size distribution: %d\n" , sd->sdCase );
01520 ShowMe();
01521 cdEXIT(EXIT_FAILURE);
01522 }
01523
01524 if( *error < 2 )
01525 ASSERT( *cs_abs > 0. && *cs_sct > 0. );
01526 if( *error < 1 )
01527 ASSERT( fabs(*cosb) <= 1.+10.*DBL_EPSILON );
01528 return;
01529 }
01530
01531
01532
01533
01534 STATIC void mie_cs(double wavlen,
01535 const sd_data *sd,
01536 const grain_data *gd,
01537 double *cs_abs,
01538 double *cs_sct,
01539 double *cosb,
01540 int *error)
01541 {
01542 bool lgOutOfBounds;
01543 long int iflag,
01544 ind;
01545 double area,
01546 aqabs,
01547 aqext,
01548 aqphas,
01549 beta,
01550 ctbrqs = -DBL_MAX,
01551 delta,
01552 frac,
01553 nim,
01554 nre,
01555 nr1,
01556 qback,
01557 qext = -DBL_MAX,
01558 qphase,
01559 qscatt = -DBL_MAX,
01560 x,
01561 xistar;
01562
01563 DEBUG_ENTRY( "mie_cs()" );
01564
01565
01566 ASSERT( wavlen > 0. );
01567 ASSERT( sd->cSize > 0. );
01568 ASSERT( gd->ndata[gd->cAxis] > 1 && (long)gd->wavlen[gd->cAxis].size() == gd->ndata[gd->cAxis] );
01569
01570
01571
01572
01573 find_arr(wavlen,gd->wavlen[gd->cAxis],gd->ndata[gd->cAxis],&ind,&lgOutOfBounds);
01574
01575 if( lgOutOfBounds )
01576 {
01577 *error = 3;
01578 *cs_abs = -1.;
01579 *cs_sct = -1.;
01580 *cosb = -2.;
01581 return;
01582 }
01583
01584 frac = (wavlen-gd->wavlen[gd->cAxis][ind])/(gd->wavlen[gd->cAxis][ind+1]-gd->wavlen[gd->cAxis][ind]);
01585 ASSERT( frac > 0.-10.*DBL_EPSILON && frac < 1.+10.*DBL_EPSILON );
01586 nre = (1.-frac)*gd->n[gd->cAxis][ind].real() + frac*gd->n[gd->cAxis][ind+1].real();
01587 ASSERT( nre > 0. );
01588 nim = (1.-frac)*gd->n[gd->cAxis][ind].imag() + frac*gd->n[gd->cAxis][ind+1].imag();
01589 ASSERT( nim > 0. );
01590 nr1 = (1.-frac)*gd->nr1[gd->cAxis][ind] + frac*gd->nr1[gd->cAxis][ind+1];
01591 ASSERT( fabs(nre-1.-nr1) < 10.*max(nre,1.)*DBL_EPSILON );
01592
01593
01594 area = PI*pow2(sd->cSize)*1.e-8;
01595
01596 x = sd->cSize/wavlen*2.*PI;
01597
01598
01599
01600
01601 sinpar(nre,nim,x,&qext,&qphase,&qscatt,&ctbrqs,&qback,&iflag);
01602
01603
01604
01605
01606 if( iflag == 0 )
01607 {
01608 *error = 0;
01609 *cs_abs = area*(qext - qscatt);
01610 *cs_sct = area*qscatt;
01611 *cosb = ctbrqs/qscatt;
01612 }
01613 else
01614 {
01615
01616 if( x >= 100. && sqrt(nr1*nr1+nim*nim) <= 0.001 )
01617 {
01618 delta = -nr1;
01619 beta = nim;
01620
01621 anomal(x,&aqext,&aqabs,&aqphas,&xistar,delta,beta);
01622
01623
01624 *error = 1;
01625 *cs_abs = area*aqabs;
01626 *cs_sct = area*(aqext - aqabs);
01627 *cosb = -2.;
01628 }
01629
01630 else
01631 {
01632 *error = 2;
01633 *cs_abs = -1.;
01634 *cs_sct = -1.;
01635 *cosb = -2.;
01636 }
01637 }
01638 if( *error < 2 )
01639 {
01640 if( *cs_abs <= 0. || *cs_sct <= 0. )
01641 {
01642 fprintf( ioQQQ, " illegal opacity found: wavl=%.4e micron," , wavlen );
01643 fprintf( ioQQQ, " abs_cs=%.2e, sct_cs=%.2e\n" , *cs_abs , *cs_sct );
01644 fprintf( ioQQQ, " please check refractive index file...\n" );
01645 cdEXIT(EXIT_FAILURE);
01646 }
01647 }
01648 if( *error < 1 )
01649 {
01650 if( fabs(*cosb) > 1.+10.*DBL_EPSILON )
01651 {
01652 fprintf( ioQQQ, " illegal asymmetry parameter found: wavl=%.4e micron," , wavlen );
01653 fprintf( ioQQQ, " cosb=%.2e\n" , *cosb );
01654 fprintf( ioQQQ, " please check refractive index file...\n" );
01655 cdEXIT(EXIT_FAILURE);
01656 }
01657 }
01658
01659 return;
01660 }
01661
01662
01663
01664
01665 STATIC void car1_fun(double wavl,
01666 const sd_data *sd,
01667 const grain_data gdArr[],
01668 double *cs_abs,
01669 double *cs_sct,
01670 double *cosb,
01671 int *error)
01672 {
01673 DEBUG_ENTRY( "car1_fun()" );
01674
01675
01676
01677
01678 const double q_gra = 0.00;
01679 const double a_xi = 50.e-4;
01680 pah1_fun(wavl,sd,&gdArr[0],cs_abs,cs_sct,cosb,error);
01681 if( sd->cSize > a_xi )
01682 {
01683 double xi_PAH = (1.-q_gra)*min(1.,pow3(a_xi/sd->cSize));
01684 double cs_abs2;
01685
01686
01687 mie_cs(wavl,sd,&gdArr[1],&cs_abs2,cs_sct,cosb,error);
01688 *cs_abs = xi_PAH*(*cs_abs) + (1.-xi_PAH)*cs_abs2;
01689 }
01690 return;
01691 }
01692
01693
01694
01695
01696
01697
01698
01699 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 };
01700 static const double pah1_wlBand[7] = { 3.3, 6.18, 7.7, 8.6, 11.3, 12.0, 13.3 };
01701 static const double pah1_width[7] = { 0.024, 0.102, 0.24, 0.168, 0.086, 0.174, 0.174 };
01702
01703 STATIC void pah1_fun(double wavl,
01704 const sd_data *sd,
01705 const grain_data *gd,
01706 double *cs_abs,
01707 double *cs_sct,
01708 double *cosb,
01709 int *error)
01710 {
01711 long int j;
01712 double cval1,
01713 pah1_fun_v,
01714 term,
01715 term1,
01716 term2,
01717 term3,
01718 x;
01719
01720 const double p1 = 4.0e-18;
01721 const double p2 = 1.1e-18;
01722 const double p3 = 3.3e-21;
01723 const double p4 = 6.0e-21;
01724 const double p5 = 2.4e-21;
01725 const double wl1a = 5.0;
01726 const double wl1b = 7.0;
01727 const double wl1c = 9.0;
01728 const double wl1d = 9.5;
01729 const double wl2a = 11.0;
01730 const double delwl2 = 0.3;
01731
01732 const double wl2b = wl2a + delwl2;
01733 const double wl2c = 15.0;
01734 const double wl3a = 3.2;
01735 const double wl3b = 3.57;
01736 const double wl3m = (wl3a+wl3b)/2.;
01737 const double wl3sig = 0.1476;
01738 const double x1 = 4.0;
01739 const double x2 = 5.9;
01740 const double x2a = RYD_INF/1.e4;
01741 const double x3 = 0.1;
01742
01743
01744 double vol = 4.*PI/3.*pow3(sd->cSize)*1.e-12;
01745
01746 double xnc = floor(vol*gd->rho/(ATOMIC_MASS_UNIT*dense.AtomicWeight[ipCARBON]));
01747
01748
01749 double xnh = floor(sqrt(6.*xnc));
01750
01751 double xnhoc = xnh/xnc;
01752
01753 double ftoc3p3 = 100.;
01754
01755 double csVal1 = 0.;
01756 double csVal2 = 0.;
01757
01758 DEBUG_ENTRY( "pah1_fun()" );
01759
01760 x = 1./wavl;
01761
01762 if( x >= x2a )
01763 {
01764
01765 double anu_ev = x/x2a*EVRYD;
01766
01767
01768 t_ADfA::Inst().set_version( PHFIT95 );
01769
01770 term1 = t_ADfA::Inst().phfit(ipHYDROGEN+1,1,1,anu_ev);
01771 term2 = 0.;
01772 for( j=1; j <= 3; ++j )
01773 term2 += t_ADfA::Inst().phfit(ipCARBON+1,6,j,anu_ev);
01774
01775 csVal1 = (xnh*term1 + xnc*term2)*1.e-18;
01776 }
01777
01778 if( x <= 2.*x2a )
01779 {
01780 cval1 = log(sqrt(xnc)*ftoc3p3/1.2328)/12.2;
01781
01782 term = pow2(min(x,x1))*(3.*x1 - 2.*min(x,x1))/pow3(x1);
01783
01784 term1 = (0.1*x + 0.41)*pow2(max(x-x2,0.));
01785
01786
01787
01788
01789
01790
01791
01792
01793
01794
01795
01796
01797
01798
01799
01800
01801
01802
01803
01804
01805
01806
01807
01808 term2 = exp(cval1*(1. - (x1/min(x,x1))));
01809
01810 term3 = p3*exp(-pow2(x/x3))*min(x,x3)/x3;
01811
01812 csVal2 = xnc*((p1*term + p2*term1)*term2 + term3);
01813 }
01814
01815 if( x2a <= x && x <= 2.*x2a )
01816 {
01817
01818 double frac = pow2(2.-x/x2a);
01819 pah1_fun_v = exp((1.-frac)*log(csVal1) + frac*log(csVal2));
01820 }
01821 else
01822 {
01823
01824 pah1_fun_v = csVal1 + csVal2;
01825 }
01826
01827
01828
01829 if( wl1a <= wavl && wl1d >= wavl )
01830 {
01831 if( wavl < wl1b )
01832 term = p4*(wavl - wl1a)/(wl1b - wl1a);
01833 else
01834 term = ( wavl > wl1c ) ? p4*(wl1d - wavl)/(wl1d - wl1c) : p4;
01835 pah1_fun_v += term*xnc;
01836 }
01837 if( wl2a <= wavl && wl2c >= wavl )
01838 {
01839 term = ( wavl < wl2b ) ? p5*((wavl - wl2a)/delwl2) : p5*sqrt(1.-pow2((wavl-wl2a)/(wl2c-wl2a)));
01840 pah1_fun_v += term*xnc;
01841 }
01842 if( wl3m-10.*wl3sig <= wavl && wavl <= wl3m+10.*wl3sig )
01843 {
01844
01845 term = 1.1*pah1_strength[0]*exp(-0.5*pow2((wavl-wl3m)/wl3sig));
01846 pah1_fun_v += term*xnh;
01847 }
01848
01849
01850 for( j=0; j < 7; j++ )
01851 {
01852 term1 = (wavl - pah1_wlBand[j])/pah1_width[j];
01853 term = 0.;
01854 if( j == 2 )
01855 {
01856
01857
01858
01859
01860 if( term1 >= -1. && term1 < -0.5 )
01861 {
01862 term = pah1_strength[j]/(3.*pah1_width[j]);
01863 term *= 2. + 2.*term1;
01864 }
01865 if( term1 >= -0.5 && term1 <= 1.5 )
01866 term = pah1_strength[j]/(3.*pah1_width[j]);
01867 if( term1 > 1.5 && term1 <= 3. )
01868 {
01869 term = pah1_strength[j]/(3.*pah1_width[j]);
01870 term *= 2. - term1*2./3.;
01871 }
01872 }
01873 else
01874 {
01875
01876
01877
01878
01879 if( term1 >= -2. && term1 < -1. )
01880 {
01881 term = pah1_strength[j]/(3.*pah1_width[j]);
01882 term *= 2. + term1;
01883 }
01884 if( term1 >= -1. && term1 <= 1. )
01885 term = pah1_strength[j]/(3.*pah1_width[j]);
01886 if( term1 > 1. && term1 <= 2. )
01887 {
01888 term = pah1_strength[j]/(3.*pah1_width[j]);
01889 term *= 2. - term1;
01890 }
01891 }
01892 if( j == 0 || j > 2 )
01893 term *= xnhoc;
01894 pah1_fun_v += term*xnc;
01895 }
01896
01897 *cs_abs = pah1_fun_v;
01898
01899
01900 *cs_sct = 0.1*pah1_fun_v;
01901 *cosb = 0.;
01902 *error = 0;
01903
01904 return;
01905 }
01906
01907
01908
01909
01910 STATIC void ld01_fun(double wavl,
01911 const sd_data *sd,
01912 const grain_data gdArr[],
01913 double *cs_abs,
01914 double *cs_sct,
01915 double *cosb,
01916 int *error)
01917 {
01918 DEBUG_ENTRY( "ld01_fun()" );
01919
01920
01921
01922
01923
01924 const double q_gra = 0.01;
01925 const double a_xi = 50.e-4;
01926 double xi_PAH, cs_abs1, cs_abs2;
01927 if( wavl >= 1./17.25 )
01928 {
01929 pah2_fun(wavl,sd,&gdArr[0],&cs_abs1,cs_sct,cosb,error);
01930 xi_PAH = (1.-q_gra)*min(1.,pow3(a_xi/sd->cSize));
01931 }
01932 else
01933 {
01934 cs_abs1 = 0.;
01935 xi_PAH = 0.;
01936 }
01937
01938
01939 mie_cs(wavl,sd,&gdArr[1],&cs_abs2,cs_sct,cosb,error);
01940 *cs_abs = xi_PAH*cs_abs1 + (1.-xi_PAH)*cs_abs2;
01941 return;
01942 }
01943
01944
01945
01946 static const double pah2_wavl[14] = { 0.0722, 0.2175, 3.3, 6.2, 7.7, 8.6, 11.3,
01947 11.9, 12.7, 16.4, 18.3, 21.2, 23.1, 26.0 };
01948 static const double pah2_width[14] = { 0.195, 0.217, 0.012, 0.032, 0.091, 0.047, 0.018,
01949 0.025, 0.024, 0.010, 0.036, 0.038, 0.046, 0.69 };
01950 static 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,
01951 7.27e-19, 1.67e-18, 5.52e-20, 6.04e-20, 1.08e-19, 2.78e-20, 1.52e-19 };
01952 static 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,
01953 6.14e-19, 1.49e-18, 5.52e-20, 6.04e-20, 1.08e-19, 2.78e-20, 1.52e-19 };
01954
01955
01956
01957 STATIC void pah2_fun(double wavl,
01958 const sd_data *sd,
01959 const grain_data *gd,
01960 double *cs_abs,
01961 double *cs_sct,
01962 double *cosb,
01963 int *error)
01964 {
01965 DEBUG_ENTRY( "pah2_fun()" );
01966
01967
01968
01969 const double E62 = 3.;
01970 const double E77 = 2.;
01971 const double E86 = 2.;
01972
01973
01974 double vol = 4.*PI/3.*pow3(sd->cSize)*1.e-12;
01975
01976 double xnc = vol*gd->rho/(ATOMIC_MASS_UNIT*dense.AtomicWeight[ipCARBON]);
01977
01978 double xnhoc;
01979
01980 if( xnc <= 25. )
01981 xnhoc = 0.5;
01982 else if( xnc <= 100. )
01983 xnhoc = 2.5/sqrt(xnc);
01984 else
01985 xnhoc = 0.25;
01986
01987 double x = 1./wavl;
01988
01989 double pah2_fun_v = 0.;
01990 if( x < 3.3 )
01991 {
01992 double M = ( xnc <= 40. ) ? 0.3*xnc : 0.4*xnc;
01993
01994 double cutoff;
01995 if( gd->charge == 0 )
01996 cutoff = 1./(3.804/sqrt(M) + 1.052);
01997 else
01998 cutoff = 1./(2.282/sqrt(M) + 0.889);
01999 double y = cutoff/wavl;
02000
02001 double cutoff_fun = atan(1.e3*pow3(y-1.)/y)/PI + 0.5;
02002
02003 pah2_fun_v = 34.58*pow( 10., -18.-3.431/x )*cutoff_fun;
02004 for( int j=2; j < 14; ++j )
02005 {
02006 double strength = ( gd->charge == 0 ) ? pah2n_strength[j] : pah2c_strength[j];
02007 if( j == 2 )
02008 strength *= xnhoc;
02009 else if( j == 3 )
02010 strength *= E62;
02011 else if( j == 4 )
02012 strength *= E77;
02013 else if( j == 5 )
02014 strength *= E86*xnhoc;
02015 else if( j == 6 || j == 7 || j == 8 )
02016 strength *= xnhoc/3.;
02017 pah2_fun_v += Drude( wavl, pah2_wavl[j], pah2_width[j], strength );
02018 }
02019 }
02020 else if( x < 5.9 )
02021 {
02022
02023 pah2_fun_v = Drude( wavl, pah2_wavl[1], pah2_width[1], pah2n_strength[1] );
02024 pah2_fun_v += (1.8687 + 0.1905*x)*1.e-18;
02025 }
02026 else if( x < 7.7 )
02027 {
02028
02029 pah2_fun_v = Drude( wavl, pah2_wavl[1], pah2_width[1], pah2n_strength[1] );
02030 double y = x - 5.9;
02031 pah2_fun_v += (1.8687 + 0.1905*x + pow2(y)*(0.4175 + 0.04370*y))*1.e-18;
02032 }
02033 else if( x < 10. )
02034 {
02035
02036 pah2_fun_v = (((-0.1057*x + 2.950)*x - 24.367)*x + 66.302)*1.e-18;
02037 }
02038 else if( x < 15. )
02039 {
02040
02041 pah2_fun_v = Drude( wavl, pah2_wavl[0], pah2_width[0], pah2n_strength[0] );
02042 pah2_fun_v += (-3. + 1.35*x)*1.e-18;
02043 }
02044 else if( x < 17.26 )
02045 {
02046
02047 pah2_fun_v = (126.0 - 6.4943*x)*1.e-18;
02048 }
02049 else
02050
02051
02052 TotalInsanity();
02053
02054
02055 pah2_fun_v *= xnc;
02056
02057 *cs_abs = pah2_fun_v;
02058
02059 *cs_sct = 0.1*pah2_fun_v;
02060 *cosb = 0.;
02061 *error = 0;
02062
02063 return;
02064 }
02065
02066
02067
02068 inline double Drude(double lambda,
02069 double lambdac,
02070 double gamma,
02071 double sigma)
02072 {
02073 double x = lambda/lambdac;
02074
02075 return 2.e-4/PI*gamma*lambdac*sigma/(pow2(x - 1./x) + pow2(gamma));
02076 }
02077
02078
02079 STATIC void tbl_fun(double wavl,
02080 const sd_data *sd,
02081 const grain_data *gd,
02082 double *cs_abs,
02083 double *cs_sct,
02084 double *cosb,
02085 int *error)
02086 {
02087 bool lgOutOfBounds;
02088 long int ind;
02089 double anu = WAVNRYD/wavl*1.e4;
02090
02091 DEBUG_ENTRY( "tbl_fun()" );
02092
02093
02094 if( sd == NULL )
02095 TotalInsanity();
02096
02099 find_arr(anu,gd->opcAnu,gd->nOpcData,&ind,&lgOutOfBounds);
02100 if( !lgOutOfBounds )
02101 {
02102 double a1g;
02103 double frac = log(anu/gd->opcAnu[ind])/log(gd->opcAnu[ind+1]/gd->opcAnu[ind]);
02104
02105 *cs_abs = exp((1.-frac)*log(gd->opcData[0][ind])+frac*log(gd->opcData[0][ind+1]));
02106 ASSERT( *cs_abs > 0. );
02107 if( gd->nOpcCols > 1 )
02108 *cs_sct = exp((1.-frac)*log(gd->opcData[1][ind])+frac*log(gd->opcData[1][ind+1]));
02109 else
02110 *cs_sct = 0.1*(*cs_abs);
02111 ASSERT( *cs_sct > 0. );
02112 if( gd->nOpcCols > 2 )
02113 a1g = exp((1.-frac)*log(gd->opcData[2][ind])+frac*log(gd->opcData[2][ind+1]));
02114 else
02115 a1g = 1.;
02116 ASSERT( a1g > 0. );
02117 *cosb = 1. - a1g;
02118 *error = 0;
02119 }
02120 else
02121 {
02122 *cs_abs = -1.;
02123 *cs_sct = -1.;
02124 *cosb = -2.;
02125 *error = 3;
02126 }
02127 return;
02128 }
02129
02130
02131 STATIC double size_distr(double size,
02132 const sd_data *sd)
02133 {
02134 bool lgOutOfBounds;
02135 long ind;
02136 double frac,
02137 res,
02138 x;
02139
02140 DEBUG_ENTRY( "size_distr()" );
02141
02142 if( size >= sd->lim[ipBLo] && size <= sd->lim[ipBHi] )
02143 switch( sd->sdCase )
02144 {
02145 case SD_SINGLE_SIZE:
02146 case SD_NR_CARBON:
02147 res = 1.;
02148 break;
02149 case SD_POWERLAW:
02150
02151 case SD_EXP_CUTOFF1:
02152 case SD_EXP_CUTOFF2:
02153 case SD_EXP_CUTOFF3:
02154
02155
02156 res = pow(size,sd->a[ipExp]);
02157 if( sd->a[ipBeta] < 0. )
02158 res /= (1. - sd->a[ipBeta]*size);
02159 else if( sd->a[ipBeta] > 0. )
02160 res *= (1. + sd->a[ipBeta]*size);
02161 if( size < sd->a[ipBLo] && sd->a[ipSLo] > 0. )
02162 res *= exp(-powi((sd->a[ipBLo]-size)/sd->a[ipSLo],nint(sd->a[ipAlpha])));
02163 if( size > sd->a[ipBHi] && sd->a[ipSHi] > 0. )
02164 res *= exp(-powi((size-sd->a[ipBHi])/sd->a[ipSHi],nint(sd->a[ipAlpha])));
02165 break;
02166 case SD_LOG_NORMAL:
02167 x = log(size/sd->a[ipGCen])/sd->a[ipGSig];
02168 res = exp(-0.5*pow2(x))/size;
02169 break;
02170 case SD_LIN_NORMAL:
02171 x = (size-sd->a[ipGCen])/sd->a[ipGSig];
02172 res = exp(-0.5*pow2(x))/size;
02173 break;
02174 case SD_TABLE:
02175 find_arr(log(size),sd->ln_a,sd->npts,&ind,&lgOutOfBounds);
02176 if( lgOutOfBounds )
02177 {
02178 fprintf( ioQQQ, " size distribution table has insufficient range\n" );
02179 fprintf( ioQQQ, " requested size: %.5f table range %.5f - %.5f\n",
02180 size, exp(sd->ln_a[0]), exp(sd->ln_a[sd->npts-1]) );
02181 cdEXIT(EXIT_FAILURE);
02182 }
02183 frac = (log(size)-sd->ln_a[ind])/(sd->ln_a[ind+1]-sd->ln_a[ind]);
02184 ASSERT( frac > 0.-10.*DBL_EPSILON && frac < 1.+10.*DBL_EPSILON );
02185 res = (1.-frac)*sd->ln_a4dNda[ind] + frac*sd->ln_a4dNda[ind+1];
02186
02187 res = exp(res)/POW4(size);
02188 break;
02189 case SD_ILLEGAL:
02190 default:
02191 fprintf( ioQQQ, " insane case for grain size distribution: %d\n" , sd->sdCase );
02192 ShowMe();
02193 cdEXIT(EXIT_FAILURE);
02194 }
02195 else
02196 res = 0.;
02197 return res;
02198 }
02199
02200
02201
02202
02203
02204
02205 STATIC double search_limit(double ref,
02206 double step,
02207 double rel_cutoff,
02208
02209 sd_data sd)
02210 {
02211 long i;
02212 double f1,
02213 f2,
02214 fmid,
02215 renorm,
02216 x1,
02217 x2 = DBL_MAX,
02218 xmid = DBL_MAX;
02219
02220
02221 const double TOLER = 1.e-6;
02222
02223 DEBUG_ENTRY( "search_limit()" );
02224
02225
02226 ASSERT( rel_cutoff > 0. && rel_cutoff < 1. );
02227
02228 if( step == 0. )
02229 {
02230 return ref;
02231 }
02232
02233
02234
02235
02236 sd.lim[ipBLo] = 0.;
02237 sd.lim[ipBHi] = DBL_MAX;
02238
02239 x1 = ref;
02240
02241 f1 = -log(rel_cutoff);
02242 renorm = f1 - log(POW4(x1)*size_distr(x1,&sd));
02243
02244
02245 f2 = 1.;
02246 for( i=0; i < 20 && f2 > 0.; ++i )
02247 {
02248 x2 = max(ref + step,SMALLEST_GRAIN);
02249 f2 = log(POW4(x2)*size_distr(x2,&sd)) + renorm;
02250 if( f2 >= 0. )
02251 {
02252 x1 = x2;
02253 f1 = f2;
02254 }
02255 step *= 2.;
02256 }
02257 if( f2 > 0. )
02258 {
02259 fprintf( ioQQQ, " Could not bracket solution\n" );
02260 cdEXIT(EXIT_FAILURE);
02261 }
02262
02263
02264 while( 2.*fabs(x1-x2)/(x1+x2) > TOLER )
02265 {
02266 xmid = (x1+x2)/2.;
02267 fmid = log(POW4(xmid)*size_distr(xmid,&sd)) + renorm;
02268
02269 if( fmid == 0. )
02270 break;
02271
02272 if( f1*fmid > 0. )
02273 {
02274 x1 = xmid;
02275 f1 = fmid;
02276 }
02277 else
02278 {
02279 x2 = xmid;
02280 f2 = fmid;
02281 }
02282 }
02283 return (x1+x2)/2.;
02284 }
02285
02286
02287 STATIC void mie_calc_ial( const grain_data *gd,
02288 long int n,
02289 vector<double>& invlen,
02290 const char *chString,
02291 bool *lgWarning)
02292 {
02293 bool lgErrorOccurred=true,
02294 lgOutOfBounds;
02295 long int i,
02296 ind,
02297 j;
02298 double frac,
02299 InvDep,
02300 nim,
02301 wavlen;
02302
02303 DEBUG_ENTRY( "mie_calc_ial()" );
02304
02305
02306 ASSERT( gd->rfiType == RFI_TABLE );
02307
02308 vector<int> ErrorIndex( rfield.nupper );
02309
02310 for( i=0; i < n; i++ )
02311 {
02312 wavlen = WAVNRYD/rfield.anu[i]*1.e4;
02313
02314 ErrorIndex[i] = 0;
02315 lgErrorOccurred = false;
02316 invlen[i] = 0.;
02317
02318 for( j=0; j < gd->nAxes; j++ )
02319 {
02320
02321 find_arr(wavlen,gd->wavlen[j],gd->ndata[j],&ind,&lgOutOfBounds);
02322 if( lgOutOfBounds )
02323 {
02324 ErrorIndex[i] = 3;
02325 lgErrorOccurred = true;
02326 invlen[i] = 0.;
02327 break;
02328 }
02329 frac = (wavlen-gd->wavlen[j][ind])/(gd->wavlen[j][ind+1]-gd->wavlen[j][ind]);
02330 nim = (1.-frac)*gd->n[j][ind].imag() + frac*gd->n[j][ind+1].imag();
02331
02332
02333 InvDep = PI4*nim/wavlen*1.e4;
02334 ASSERT( InvDep > 0. );
02335
02336 invlen[i] += InvDep*gd->wt[j];
02337 }
02338 }
02339
02340 if( lgErrorOccurred )
02341 {
02342 mie_repair(chString,n,3,3,rfield.anu,&invlen[0],ErrorIndex,false,lgWarning);
02343 }
02344
02345 return;
02346 }
02347
02348
02349
02350 const int NPTS_DERIV = 8;
02351 const int NPTS_COMB = NPTS_DERIV*(NPTS_DERIV-1)/2;
02352
02353
02354 STATIC void mie_repair( const char *chString,
02355 long int n,
02356 int val,
02357 int del,
02358 const realnum anu[],
02359 double data[],
02360 vector<int>& ErrorIndex,
02361 bool lgRound,
02362 bool *lgWarning)
02363 {
02364 bool lgExtrapolate,
02365 lgVerbose;
02366 long int i1,
02367 i2,
02368 ind1,
02369 ind2,
02370 j;
02371 double dx,
02372 sgn,
02373 slp1,
02374 xlg1,
02375 xlg2,
02376 y1lg1,
02377 y1lg2;
02378
02379
02380 const long BIG_INTERPOLATION = 10;
02381
02382 DEBUG_ENTRY( "mie_repair()" );
02383
02384 lgVerbose = ( chString[0] != '\0' );
02385
02386 for( ind1=0; ind1 < n; )
02387 {
02388 if( ErrorIndex[ind1] == val )
02389 {
02390
02391 ind2 = ind1;
02392 while( ind2 < n && ErrorIndex[ind2] == val )
02393 ind2++;
02394
02395 if( lgVerbose )
02396 fprintf( ioQQQ, " %s", chString );
02397
02398 if( ind1 == 0 )
02399 {
02400
02401 i1 = ind2;
02402 i2 = ind2+NPTS_DERIV-1;
02403 lgExtrapolate = true;
02404 sgn = +1.;
02405 if( lgVerbose )
02406 {
02407 fprintf( ioQQQ, " extrapolated below %.4e Ryd\n",anu[i1] );
02408 }
02409 }
02410 else if( ind2 == n )
02411 {
02412
02413 i1 = ind1-NPTS_DERIV;
02414 i2 = ind1-1;
02415 lgExtrapolate = true;
02416 sgn = -1.;
02417 if( lgVerbose )
02418 {
02419 fprintf( ioQQQ, " extrapolated above %.4e Ryd\n",anu[i2] );
02420 }
02421 }
02422 else
02423 {
02424
02425 i1 = ind1-1;
02426 i2 = ind2;
02427 lgExtrapolate = false;
02428 sgn = 0.;
02429 if( lgVerbose )
02430 {
02431 fprintf( ioQQQ, " interpolated between %.4e and %.4e Ryd\n",
02432 anu[i1],anu[i2] );
02433 }
02434 if( i2-i1-1 > BIG_INTERPOLATION )
02435 {
02436 if( lgVerbose )
02437 {
02438 fprintf( ioQQQ, " ***Warning: extensive interpolation used\n" );
02439 }
02440 *lgWarning = true;
02441 }
02442 }
02443
02444 if( i1 < 0 || i2 >= n )
02445 {
02446 fprintf( ioQQQ, " Insufficient data for extrapolation\n" );
02447 cdEXIT(EXIT_FAILURE);
02448 }
02449
02450 xlg1 = log(anu[i1]);
02451 y1lg1 = log(data[i1]);
02452
02453 if( lgExtrapolate )
02454 slp1 = mie_find_slope(anu,data,ErrorIndex,i1,i2,val,lgVerbose,lgWarning);
02455 else
02456 {
02457 xlg2 = log(anu[i2]);
02458 y1lg2 = log(data[i2]);
02459 slp1 = (y1lg2-y1lg1)/(xlg2-xlg1);
02460 }
02461 if( lgRound && lgExtrapolate && sgn > 0. )
02462 {
02463
02464
02465
02466 slp1 = max(slp1,0.);
02467 }
02468
02469 else if( lgExtrapolate && sgn*slp1 < 0. )
02470 {
02471 fprintf( ioQQQ, " Illegal value for slope in extrapolation %.6e\n", slp1 );
02472 cdEXIT(EXIT_FAILURE);
02473 }
02474
02475 for( j=ind1; j < ind2; j++ )
02476 {
02477 dx = log(anu[j]) - xlg1;
02478 data[j] = exp(y1lg1 + dx*slp1);
02479 ErrorIndex[j] -= del;
02480 }
02481
02482 ind1 = ind2;
02483 }
02484 else
02485 {
02486 ind1++;
02487 }
02488 }
02489
02490 for( j=0; j < n; j++ )
02491 {
02492 if( ErrorIndex[j] > val-del )
02493 {
02494 fprintf( ioQQQ, " Internal error in mie_repair, index=%ld, val=%d\n",j,ErrorIndex[j] );
02495 ShowMe();
02496 cdEXIT(EXIT_FAILURE);
02497 }
02498 }
02499 return;
02500 }
02501
02502
02503 STATIC double mie_find_slope( const realnum anu[],
02504 const double data[],
02505 const vector<int>& ErrorIndex,
02506 long i1,
02507 long i2,
02508 int val,
02509 bool lgVerbose,
02510 bool *lgWarning)
02511 {
02512 long i,
02513 j,
02514 k;
02515 double s1,
02516 s2,
02517 slope,
02518 slp1[NPTS_COMB],
02519 stdev;
02520
02521
02522
02523 const double LARGE_STDEV = 0.2;
02524
02525 DEBUG_ENTRY( "mie_find_slope()" );
02526
02527
02528 ASSERT( i2-i1 == NPTS_DERIV-1 );
02529 for( i=i1; i <= i2; i++ )
02530 {
02531 ASSERT( ErrorIndex[i] < val );
02532 ASSERT( anu[i] > 0. && data[i] > 0. );
02533 }
02534
02535
02536 for( i=0; i < NPTS_COMB; i++ )
02537 {
02538 slp1[i] = -DBL_MAX;
02539 }
02540
02541 k = 0;
02542
02543
02544 for( i=i1; i < i2; i++ )
02545 {
02546 for( j=i+1; j <= i2; j++ )
02547 {
02548 slp1[k++] = log(data[j]/data[i])/log(anu[j]/anu[i]);
02549 }
02550 }
02551
02552 for( i=0; i <= NPTS_COMB/2; i++ )
02553 {
02554 for( j=i+1; j < NPTS_COMB; j++ )
02555 {
02556 if( slp1[i] > slp1[j] )
02557 {
02558 double xxx = slp1[i];
02559 slp1[i] = slp1[j];
02560 slp1[j] = xxx;
02561 }
02562 }
02563 }
02564
02565 slope = ( NPTS_COMB%2 == 1 ) ? slp1[NPTS_COMB/2] : (slp1[NPTS_COMB/2-1] + slp1[NPTS_COMB/2])/2.;
02566
02567
02568 s1 = s2 = 0.;
02569 for( i=0; i < NPTS_COMB; i++ )
02570 {
02571 s1 += slp1[i];
02572 s2 += pow2(slp1[i]);
02573 }
02574
02575 stdev = max(s2/(double)NPTS_COMB - pow2(s1/(double)NPTS_COMB),0.);
02576 stdev = sqrt(stdev);
02577
02578
02579 if( stdev > LARGE_STDEV )
02580 {
02581 if( lgVerbose )
02582 {
02583 fprintf( ioQQQ, " ***Warning: slope for extrapolation may be unreliable\n" );
02584 }
02585 *lgWarning = true;
02586 }
02587 return slope;
02588 }
02589
02590
02591
02592 STATIC void mie_read_rfi( const char *chFile,
02593 grain_data *gd)
02594 {
02595 bool lgLogData=false;
02596 long int dl,
02597 help,
02598 i,
02599 nelem,
02600 j,
02601 nridf,
02602 sgn = 0;
02603 double eps1,
02604 eps2,
02605 LargestLog,
02606 molw,
02607 nAtoms,
02608 nr,
02609 ni,
02610 tmp1,
02611 tmp2,
02612 total = 0.;
02613 char chLine[FILENAME_PATH_LENGTH_2],
02614 chWord[FILENAME_PATH_LENGTH_2];
02615 FILE *io2;
02616
02617 DEBUG_ENTRY( "mie_read_rfi()" );
02618
02619 io2 = open_data( chFile, "r", AS_LOCAL_ONLY );
02620
02621 dl = 0;
02622
02623
02624 mie_next_data(chFile,io2,chLine,&dl);
02625 mie_read_long(chFile,chLine,&gd->magic,true,dl);
02626 if( gd->magic != MAGIC_RFI )
02627 {
02628 fprintf( ioQQQ, " Refractive index file %s has obsolete magic number\n",chFile );
02629 fprintf( ioQQQ, " I found %ld, but expected %ld on line #%ld\n",gd->magic,MAGIC_RFI,dl );
02630 fprintf( ioQQQ, " Please replace this file with an up to date version\n" );
02631 cdEXIT(EXIT_FAILURE);
02632 }
02633
02634
02635 mie_next_data(chFile,io2,chLine,&dl);
02636 mie_read_word(chLine,chWord,FILENAME_PATH_LENGTH_2,false);
02637 mie_read_form(chWord,gd->elmAbun,&nAtoms,&molw);
02638
02639
02640 gd->mol_weight = molw;
02641 gd->atom_weight = gd->mol_weight/nAtoms;
02642
02643
02644 mie_next_data(chFile,io2,chLine,&dl);
02645 mie_read_double(chFile,chLine,&gd->abun,true,dl);
02646
02647
02648 mie_next_data(chFile,io2,chLine,&dl);
02649 mie_read_double(chFile,chLine,&gd->depl,true,dl);
02650 if( gd->depl > 1. )
02651 {
02652 fprintf( ioQQQ, " Illegal value for default depletion in %s\n",chFile );
02653 fprintf( ioQQQ, " Line #%ld, depl=%14.6e\n",dl,gd->depl);
02654 cdEXIT(EXIT_FAILURE);
02655 }
02656
02657 for( nelem=0; nelem < LIMELM; nelem++ )
02658 gd->elmAbun[nelem] *= gd->abun*gd->depl;
02659
02660
02661 mie_next_data(chFile,io2,chLine,&dl);
02662 mie_read_double(chFile,chLine,&gd->rho,false,dl);
02663
02664
02665 mie_next_data(chFile,io2,chLine,&dl);
02666 mie_read_long(chFile,chLine,&help,true,dl);
02667 gd->matType = (mat_type)help;
02668 if( gd->matType >= MAT_TOP )
02669 {
02670 fprintf( ioQQQ, " Illegal value for material type in %s\n",chFile );
02671 fprintf( ioQQQ, " Line #%ld, type=%d\n",dl,gd->matType);
02672 cdEXIT(EXIT_FAILURE);
02673 }
02674
02675
02676 mie_next_data(chFile,io2,chLine,&dl);
02677 mie_read_double(chFile,chLine,&gd->work,true,dl);
02678
02679
02680 mie_next_data(chFile,io2,chLine,&dl);
02681 mie_read_double(chFile,chLine,&gd->bandgap,false,dl);
02682 if( gd->bandgap >= gd->work )
02683 {
02684 fprintf( ioQQQ, " Illegal value for bandgap in %s\n",chFile );
02685 fprintf( ioQQQ, " Line #%ld, bandgap=%.4e, work function=%.4e\n",dl,gd->bandgap,gd->work);
02686 fprintf( ioQQQ, " Bandgap should always be less than work function\n");
02687 cdEXIT(EXIT_FAILURE);
02688 }
02689
02690
02691 mie_next_data(chFile,io2,chLine,&dl);
02692 mie_read_double(chFile,chLine,&gd->therm_eff,true,dl);
02693 if( gd->therm_eff > 1.f )
02694 {
02695 fprintf( ioQQQ, " Illegal value for thermionic efficiency in %s\n",chFile );
02696 fprintf( ioQQQ, " Line #%ld, value=%.4e\n",dl,gd->therm_eff);
02697 fprintf( ioQQQ, " Allowed values are 0. < efficiency <= 1.\n");
02698 cdEXIT(EXIT_FAILURE);
02699 }
02700
02701
02702 mie_next_data(chFile,io2,chLine,&dl);
02703 mie_read_double(chFile,chLine,&gd->subl_temp,true,dl);
02704
02705
02706 mie_next_data(chFile,io2,chLine,&dl);
02707 mie_read_word(chLine,chWord,WORDLEN,true);
02708
02709 if( nMatch( "RFI_", chWord ) )
02710 gd->rfiType = RFI_TABLE;
02711 else if( nMatch( "OPC_", chWord ) )
02712 gd->rfiType = OPC_TABLE;
02713 else if( nMatch( "GREY", chWord ) )
02714 gd->rfiType = OPC_GREY;
02715 else if( nMatch( "PAH1", chWord ) )
02716 gd->rfiType = OPC_PAH1;
02717 else if( nMatch( "PH2N", chWord ) )
02718 gd->rfiType = OPC_PAH2N;
02719 else if( nMatch( "PH2C", chWord ) )
02720 gd->rfiType = OPC_PAH2C;
02721 else
02722 {
02723 fprintf( ioQQQ, " Illegal keyword in %s\n",chFile );
02724 fprintf( ioQQQ, " Line #%ld, value=%s\n",dl,chWord);
02725 fprintf( ioQQQ, " Allowed values are: RFI_TBL, OPC_TBL, GREY, PAH1, PH2N, PH2C\n");
02726 cdEXIT(EXIT_FAILURE);
02727 }
02728
02729 switch( gd->rfiType )
02730 {
02731 case RFI_TABLE:
02732
02733
02734 mie_next_data(chFile,io2,chLine,&dl);
02735 mie_read_long(chFile,chLine,&nridf,true,dl);
02736 if( nridf > 3 )
02737 {
02738 fprintf( ioQQQ, " Illegal data code in %s\n",chFile );
02739 fprintf( ioQQQ, " Line #%ld, data code=%ld\n",dl,nridf);
02740 cdEXIT(EXIT_FAILURE);
02741 }
02742
02743
02744
02745 mie_next_data(chFile,io2,chLine,&dl);
02746 mie_read_long(chFile,chLine,&gd->nAxes,true,dl);
02747 if( gd->nAxes > NAX )
02748 {
02749 fprintf( ioQQQ, " Illegal no. of axes in %s\n",chFile );
02750 fprintf( ioQQQ, " Line #%ld, number=%ld\n",dl,gd->nAxes);
02751 cdEXIT(EXIT_FAILURE);
02752 }
02753
02754
02755 mie_next_data(chFile,io2,chLine,&dl);
02756 switch( gd->nAxes )
02757 {
02758 case 1:
02759 mie_read_double(chFile,chLine,&gd->wt[0],true,dl);
02760 total = gd->wt[0];
02761 break;
02762 case 2:
02763 if( sscanf( chLine, "%lf %lf", &gd->wt[0], &gd->wt[1] ) != 2 )
02764 {
02765 fprintf( ioQQQ, " Syntax error in %s\n",chFile);
02766 fprintf( ioQQQ, " Line #%ld: %s\n",dl,chLine);
02767 cdEXIT(EXIT_FAILURE);
02768 }
02769 if( gd->wt[0] <= 0. || gd->wt[1] <= 0. )
02770 {
02771 fprintf( ioQQQ, " Illegal data in %s\n",chFile);
02772 fprintf( ioQQQ, " Line #%ld: %s\n",dl,chLine);
02773 cdEXIT(EXIT_FAILURE);
02774 }
02775 total = gd->wt[0] + gd->wt[1];
02776 break;
02777 case 3:
02778 if( sscanf( chLine, "%lf %lf %lf", &gd->wt[0], &gd->wt[1], &gd->wt[2] ) != 3 )
02779 {
02780 fprintf( ioQQQ, " Syntax error in %s\n",chFile);
02781 fprintf( ioQQQ, " Line #%ld: %s\n",dl,chLine);
02782 cdEXIT(EXIT_FAILURE);
02783 }
02784 if( gd->wt[0] <= 0. || gd->wt[1] <= 0. || gd->wt[2] <= 0. )
02785 {
02786 fprintf( ioQQQ, " Illegal data in %s\n",chFile);
02787 fprintf( ioQQQ, " Line #%ld: %s\n",dl,chLine);
02788 cdEXIT(EXIT_FAILURE);
02789 }
02790 total = gd->wt[0] + gd->wt[1] + gd->wt[2];
02791 break;
02792 default:
02793 fprintf( ioQQQ, " insane case for gd->nAxes: %ld\n", gd->nAxes );
02794 ShowMe();
02795 cdEXIT(EXIT_FAILURE);
02796 }
02797 for( j=0; j < gd->nAxes; j++ )
02798 {
02799 gd->wt[j] /= total;
02800
02801
02802 mie_next_data(chFile,io2,chLine,&dl);
02803 mie_read_long(chFile,chLine,&gd->ndata[j],false,dl);
02804 if( gd->ndata[j] < 2 )
02805 {
02806 fprintf( ioQQQ, " Illegal number of data points in %s\n",chFile );
02807 fprintf( ioQQQ, " Line #%ld, number=%ld\n",dl,gd->ndata[j]);
02808 cdEXIT(EXIT_FAILURE);
02809 }
02810
02811
02812 gd->wavlen[j].resize( gd->ndata[j] );
02813 gd->n[j].resize( gd->ndata[j] );
02814 gd->nr1[j].resize( gd->ndata[j] );
02815
02816 for( i=0; i < gd->ndata[j]; i++ )
02817 {
02818
02819
02820 mie_next_data(chFile,io2,chLine,&dl);
02821 if( sscanf( chLine, "%lf %lf %lf", &gd->wavlen[j][i], &nr, &ni ) != 3 )
02822 {
02823 fprintf( ioQQQ, " Syntax error in %s\n",chFile);
02824 fprintf( ioQQQ, " Line #%ld: %s\n",dl,chLine);
02825 cdEXIT(EXIT_FAILURE);
02826 }
02827 if( gd->wavlen[j][i] <= 0. )
02828 {
02829 fprintf( ioQQQ, " Illegal value for wavelength in %s\n",chFile);
02830 fprintf( ioQQQ, " Line #%ld, wavl=%14.6e\n",dl,gd->wavlen[j][i]);
02831 cdEXIT(EXIT_FAILURE);
02832 }
02833
02834
02835 if( i == 1 )
02836 {
02837 sgn = sign3(gd->wavlen[j][1]-gd->wavlen[j][0]);
02838 if( sgn == 0 )
02839 {
02840 fprintf( ioQQQ, " Illegal value for wavelength in %s\n",chFile);
02841 fprintf( ioQQQ, " Line #%ld, wavl=%14.6e\n",dl,gd->wavlen[j][i]);
02842 cdEXIT(EXIT_FAILURE);
02843 }
02844 }
02845 else if( i > 1 )
02846 {
02847 if( sign3(gd->wavlen[j][i]-gd->wavlen[j][i-1]) != sgn )
02848 {
02849 fprintf( ioQQQ, " Illegal value for wavelength in %s\n",chFile);
02850 fprintf( ioQQQ, " Line #%ld, wavl=%14.6e\n",dl,gd->wavlen[j][i]);
02851 cdEXIT(EXIT_FAILURE);
02852 }
02853 }
02854
02855
02856
02857 switch( nridf )
02858 {
02859 case 1:
02860 eps1 = nr;
02861 eps2 = ni;
02862 dftori(&nr,&ni,eps1,eps2);
02863 gd->nr1[j][i] = nr - 1.;
02864 break;
02865 case 2:
02866 gd->nr1[j][i] = nr;
02867 nr += 1.;
02868 break;
02869 case 3:
02870 gd->nr1[j][i] = nr - 1.;
02871 break;
02872 default:
02873 fprintf( ioQQQ, " insane case for nridf: %ld\n", nridf );
02874 ShowMe();
02875 cdEXIT(EXIT_FAILURE);
02876 }
02877 gd->n[j][i] = complex<double>(nr,ni);
02878
02879
02880 if( nr <= 0. || ni < 0. )
02881 {
02882 fprintf( ioQQQ, " Illegal value for refractive index in %s\n",chFile);
02883 fprintf( ioQQQ, " Line #%ld, (nr,ni)=(%14.6e,%14.6e)\n",dl,nr,ni);
02884 cdEXIT(EXIT_FAILURE);
02885 }
02886 ASSERT( fabs(nr-1.-gd->nr1[j][i]) < 10.*nr*DBL_EPSILON );
02887 }
02888 }
02889 break;
02890 case OPC_TABLE:
02891
02892
02893
02894
02895
02896
02897 mie_next_data(chFile,io2,chLine,&dl);
02898 mie_read_long(chFile,chLine,&gd->nOpcCols,true,dl);
02899 if( gd->nOpcCols > NDAT )
02900 {
02901 fprintf( ioQQQ, " Illegal no. of data columns in %s\n",chFile );
02902 fprintf( ioQQQ, " Line #%ld, number=%ld\n",dl,gd->nOpcCols);
02903 cdEXIT(EXIT_FAILURE);
02904 }
02905
02906
02907 mie_next_data(chFile,io2,chLine,&dl);
02908 mie_read_word(chLine,chWord,WORDLEN,true);
02909
02910 if( nMatch( "LINE", chWord ) )
02911 lgLogData = false;
02912 else if( nMatch( "LOG", chWord ) )
02913 lgLogData = true;
02914 else
02915 {
02916 fprintf( ioQQQ, " Keyword not recognized in %s\n",chFile );
02917 fprintf( ioQQQ, " Line #%ld, keyword=%s\n",dl,chWord);
02918 cdEXIT(EXIT_FAILURE);
02919 }
02920
02921
02922
02923 mie_next_data(chFile,io2,chLine,&dl);
02924 mie_read_long(chFile,chLine,&gd->nOpcData,false,dl);
02925 if( gd->nOpcData < 2 )
02926 {
02927 fprintf( ioQQQ, " Illegal number of data points in %s\n",chFile );
02928 fprintf( ioQQQ, " Line #%ld, number=%ld\n",dl,gd->nOpcData);
02929 cdEXIT(EXIT_FAILURE);
02930 }
02931
02932
02933 gd->opcAnu.resize(gd->nOpcData);
02934 for( j=0; j < gd->nOpcCols; j++ )
02935 gd->opcData[j].resize(gd->nOpcData);
02936
02937 tmp1 = -log10(1.01*DBL_MIN);
02938 tmp2 = log10(0.99*DBL_MAX);
02939 LargestLog = min(tmp1,tmp2);
02940
02941
02942
02943
02944
02945
02946
02947
02948
02949
02950 for( i=0; i < gd->nOpcData; i++ )
02951 {
02952 mie_next_data(chFile,io2,chLine,&dl);
02953 switch( gd->nOpcCols )
02954 {
02955 case 1:
02956 if( sscanf( chLine, "%lf %lf", &gd->opcAnu[i], &gd->opcData[0][i] ) != 2 )
02957 {
02958 fprintf( ioQQQ, " Syntax error in %s\n",chFile);
02959 fprintf( ioQQQ, " Line #%ld: %s\n",dl,chLine);
02960 cdEXIT(EXIT_FAILURE);
02961 }
02962 break;
02963 case 2:
02964 if( sscanf( chLine, "%lf %lf %lf", &gd->opcAnu[i], &gd->opcData[0][i],
02965 &gd->opcData[1][i] ) != 3 )
02966 {
02967 fprintf( ioQQQ, " Syntax error in %s\n",chFile);
02968 fprintf( ioQQQ, " Line #%ld: %s\n",dl,chLine);
02969 cdEXIT(EXIT_FAILURE);
02970 }
02971 break;
02972 case 3:
02973 if( sscanf( chLine, "%lf %lf %lf %lf", &gd->opcAnu[i], &gd->opcData[0][i],
02974 &gd->opcData[1][i], &gd->opcData[2][i] ) != 4 )
02975 {
02976 fprintf( ioQQQ, " Syntax error in %s\n",chFile);
02977 fprintf( ioQQQ, " Line #%ld: %s\n",dl,chLine);
02978 cdEXIT(EXIT_FAILURE);
02979 }
02980 break;
02981 case 4:
02982 if( sscanf( chLine, "%lf %lf %lf %lf %lf", &gd->opcAnu[i], &gd->opcData[0][i],
02983 &gd->opcData[1][i], &gd->opcData[2][i], &gd->opcData[3][i] ) != 5 )
02984 {
02985 fprintf( ioQQQ, " Syntax error in %s\n",chFile);
02986 fprintf( ioQQQ, " Line #%ld: %s\n",dl,chLine);
02987 cdEXIT(EXIT_FAILURE);
02988 }
02989 break;
02990 default:
02991 fprintf( ioQQQ, " insane case for gd->nOpcCols: %ld\n", gd->nOpcCols );
02992 ShowMe();
02993 cdEXIT(EXIT_FAILURE);
02994 }
02995 if( lgLogData )
02996 {
02997
02998 if( fabs(gd->opcAnu[i]) > LargestLog )
02999 {
03000 fprintf( ioQQQ, " Illegal value for frequency in %s\n",chFile );
03001 fprintf( ioQQQ, " Line #%ld, freq=%14.6e\n",dl,gd->opcAnu[i] );
03002 cdEXIT(EXIT_FAILURE);
03003 }
03004 gd->opcAnu[i] = pow(10.,gd->opcAnu[i]);
03005 for( j=0; j < gd->nOpcCols; j++ )
03006 {
03007 if( fabs(gd->opcData[j][i]) > LargestLog )
03008 {
03009 fprintf( ioQQQ, " Illegal data value in %s\n",chFile );
03010 fprintf( ioQQQ, " Line #%ld, value=%14.6e\n",dl,gd->opcData[j][i] );
03011 cdEXIT(EXIT_FAILURE);
03012 }
03013 gd->opcData[j][i] = pow(10.,gd->opcData[j][i]);
03014 }
03015 }
03016 if( gd->opcAnu[i] <= 0. )
03017 {
03018 fprintf( ioQQQ, " Illegal value for frequency in %s\n",chFile );
03019 fprintf( ioQQQ, " Line #%ld, freq=%14.6e\n",dl,gd->opcAnu[i] );
03020 cdEXIT(EXIT_FAILURE);
03021 }
03022 for( j=0; j < gd->nOpcCols; j++ )
03023 {
03024 if( gd->opcData[j][i] <= 0. )
03025 {
03026 fprintf( ioQQQ, " Illegal data value in %s\n",chFile );
03027 fprintf( ioQQQ, " Line #%ld, value=%14.6e\n",dl,gd->opcData[j][i] );
03028 cdEXIT(EXIT_FAILURE);
03029 }
03030 }
03031
03032
03033 if( i == 1 )
03034 {
03035 sgn = sign3(gd->opcAnu[1]-gd->opcAnu[0]);
03036 if( sgn == 0 )
03037 {
03038 double dataVal = lgLogData ? log10(gd->opcAnu[1]) : gd->opcAnu[1];
03039 fprintf( ioQQQ, " Illegal value for frequency in %s\n",chFile );
03040 fprintf( ioQQQ, " Line #%ld, freq=%14.6e\n",dl,dataVal );
03041 cdEXIT(EXIT_FAILURE);
03042 }
03043 }
03044 else if( i > 1 )
03045 {
03046 if( sign3(gd->opcAnu[i]-gd->opcAnu[i-1]) != sgn )
03047 {
03048 double dataVal = lgLogData ? log10(gd->opcAnu[i]) : gd->opcAnu[i];
03049 fprintf( ioQQQ, " Illegal value for frequency in %s\n",chFile);
03050 fprintf( ioQQQ, " Line #%ld, freq=%14.6e\n",dl,dataVal);
03051 cdEXIT(EXIT_FAILURE);
03052 }
03053 }
03054 }
03055 gd->nAxes = 1;
03056 break;
03057 case OPC_GREY:
03058 case OPC_PAH1:
03059 case OPC_PAH2N:
03060 case OPC_PAH2C:
03061
03062
03063 gd->nAxes = 1;
03064 break;
03065 default:
03066 fprintf( ioQQQ, " Insane value for gd->rfiType: %d, bailing out\n", gd->rfiType );
03067 ShowMe();
03068 cdEXIT(EXIT_FAILURE);
03069 }
03070
03071 fclose(io2);
03072 return;
03073 }
03074
03075
03076
03077 STATIC void mie_read_mix( const char *chFile,
03078 grain_data *gd)
03079 {
03080 emt_type EMTtype;
03081 long int dl,
03082 i,
03083 j,
03084 k,
03085 l,
03086 nelem,
03087 nMaterial,
03088 sumAxes;
03089 double maxIndex = DBL_MAX,
03090 minIndex = DBL_MAX,
03091 nAtoms,
03092 sum,
03093 sum2,
03094 wavHi,
03095 wavLo,
03096 wavStep;
03097 complex<double> eps_eff(-DBL_MAX,-DBL_MAX);
03098 char chLine[FILENAME_PATH_LENGTH_2],
03099 chWord[FILENAME_PATH_LENGTH_2],
03100 *str;
03101 FILE *io2;
03102
03103 DEBUG_ENTRY( "mie_read_mix()" );
03104
03105 io2 = open_data( chFile, "r", AS_LOCAL_ONLY );
03106
03107 dl = 0;
03108
03109
03110 mie_next_data(chFile,io2,chLine,&dl);
03111 mie_read_long(chFile,chLine,&gd->magic,true,dl);
03112 if( gd->magic != MAGIC_MIX )
03113 {
03114 fprintf( ioQQQ, " Mixed grain file %s has obsolete magic number\n",chFile );
03115 fprintf( ioQQQ, " I found %ld, but expected %ld on line #%ld\n",gd->magic,MAGIC_MIX,dl );
03116 fprintf( ioQQQ, " Please replace this file with an up to date version\n" );
03117 cdEXIT(EXIT_FAILURE);
03118 }
03119
03120
03121 mie_next_data(chFile,io2,chLine,&dl);
03122 mie_read_double(chFile,chLine,&gd->depl,true,dl);
03123 if( gd->depl > 1. )
03124 {
03125 fprintf( ioQQQ, " Illegal value for default depletion in %s\n",chFile );
03126 fprintf( ioQQQ, " Line #%ld, depl=%14.6e\n",dl,gd->depl);
03127 cdEXIT(EXIT_FAILURE);
03128 }
03129
03130
03131 mie_next_data(chFile,io2,chLine,&dl);
03132 mie_read_long(chFile,chLine,&nMaterial,true,dl);
03133 if( nMaterial < 2 )
03134 {
03135 fprintf( ioQQQ, " Illegal number of materials in mixed grain file %s\n",chFile );
03136 fprintf( ioQQQ, " I found %ld on line #%ld\n",nMaterial,dl );
03137 fprintf( ioQQQ, " This number should be at least 2\n" );
03138 cdEXIT(EXIT_FAILURE);
03139 }
03140
03141 vector<double> frac(nMaterial);
03142 vector<double> frac2(nMaterial);
03143 vector<grain_data> gdArr(nMaterial);
03144
03145 sum = 0.;
03146 sum2 = 0.;
03147 sumAxes = 0;
03148 for( i=0; i < nMaterial; i++ )
03149 {
03150 char chFile2[FILENAME_PATH_LENGTH_2];
03151
03152
03153
03154 mie_next_data(chFile,io2,chLine,&dl);
03155 mie_read_double(chFile,chLine,&frac[i],true,dl);
03156
03157 sum += frac[i];
03158
03159 str = strchr_s( chLine, '\"' );
03160 if( str != NULL )
03161 {
03162 strcpy( chFile2, ++str );
03163 str = strchr_s( chFile2, '\"' );
03164 if( str != NULL )
03165 *str = '\0';
03166 }
03167 if( str == NULL )
03168 {
03169 fprintf( ioQQQ, " No pair of double quotes was found on line #%ld of file %s\n",dl,chFile );
03170 fprintf( ioQQQ, " Please supply the refractive index file name between double quotes\n" );
03171 cdEXIT(EXIT_FAILURE);
03172 }
03173
03174 mie_read_rfi( chFile2, &gdArr[i] );
03175 if( gdArr[i].rfiType != RFI_TABLE )
03176 {
03177 fprintf( ioQQQ, " Input error on line #%ld of file %s\n",dl,chFile );
03178 fprintf( ioQQQ, " File %s is not of type RFI_TBL, this is illegal\n",chFile2 );
03179 cdEXIT(EXIT_FAILURE);
03180 }
03181
03182
03183 if( i == nMaterial-1 && gdArr[i].mol_weight == 0. )
03184 {
03185 fprintf( ioQQQ, " Input error on line #%ld of file %s\n",dl,chFile );
03186 fprintf( ioQQQ, " Last entry in list of materials is vacuum, this is not allowed\n" );
03187 fprintf( ioQQQ, " Please move this entry to an earlier position in the file\n" );
03188 cdEXIT(EXIT_FAILURE);
03189 }
03190
03191 frac2[i] = ( gdArr[i].mol_weight > 0. ) ? frac[i] : 0.;
03192 sum2 += frac2[i];
03193
03194 sumAxes += gdArr[i].nAxes;
03195 }
03196
03197
03198 mie_next_data(chFile,io2,chLine,&dl);
03199 mie_read_word(chLine,chWord,WORDLEN,true);
03200
03201 if( nMatch( "FA00", chWord ) )
03202 EMTtype = FARAFONOV00;
03203 else if( nMatch( "ST95", chWord ) )
03204 EMTtype = STOGNIENKO95;
03205 else if( nMatch( "BR35", chWord ) )
03206 EMTtype = BRUGGEMAN35;
03207 else
03208 {
03209 fprintf( ioQQQ, " Keyword not recognized in %s\n",chFile );
03210 fprintf( ioQQQ, " Line #%ld, keyword=%s\n",dl,chWord);
03211 cdEXIT(EXIT_FAILURE);
03212 }
03213
03214
03215 for( i=0; i < nMaterial; i++ )
03216 {
03217 frac[i] /= sum;
03218 frac2[i] /= sum2;
03219
03220 for( nelem=0; nelem < LIMELM; nelem++ )
03221 {
03222 gdArr[i].elmAbun[nelem] /= gdArr[i].abun*gdArr[i].depl;
03223 }
03224 }
03225
03226 wavLo = 0.;
03227 wavHi = DBL_MAX;
03228 gd->abun = DBL_MAX;
03229 for( nelem=0; nelem < LIMELM; nelem++ )
03230 {
03231 gd->elmAbun[nelem] = 0.;
03232 }
03233 gd->mol_weight = 0.;
03234 gd->rho = 0.;
03235 gd->work = DBL_MAX;
03236 gd->bandgap = DBL_MAX;
03237 gd->therm_eff = 0.;
03238 gd->subl_temp = DBL_MAX;
03239 gd->nAxes = 1;
03240 gd->wt[0] = 1.;
03241 gd->ndata[0] = MIX_TABLE_SIZE;
03242 gd->rfiType = RFI_TABLE;
03243
03244 for( i=0; i < nMaterial; i++ )
03245 {
03246 for( k=0; k < gdArr[i].nAxes; k++ )
03247 {
03248 double wavMin = min(gdArr[i].wavlen[k][0],gdArr[i].wavlen[k][gdArr[i].ndata[k]-1]);
03249 double wavMax = max(gdArr[i].wavlen[k][0],gdArr[i].wavlen[k][gdArr[i].ndata[k]-1]);
03250 wavLo = max(wavLo,wavMin);
03251 wavHi = min(wavHi,wavMax);
03252 }
03253 minIndex = DBL_MAX;
03254 maxIndex = 0.;
03255 for( nelem=0; nelem < LIMELM; nelem++ )
03256 {
03257 gd->elmAbun[nelem] += frac[i]*gdArr[i].elmAbun[nelem];
03258 if( gd->elmAbun[nelem] > 0. )
03259 {
03260 minIndex = min(minIndex,gd->elmAbun[nelem]);
03261 }
03262 maxIndex = max(maxIndex,gd->elmAbun[nelem]);
03263 }
03264 gd->mol_weight += frac[i]*gdArr[i].mol_weight;
03265 gd->rho += frac[i]*gdArr[i].rho;
03266
03267 if( gdArr[i].mol_weight > 0. )
03268 {
03269 gd->abun = min(gd->abun,gdArr[i].abun/frac[i]);
03270 switch( EMTtype )
03271 {
03272 case FARAFONOV00:
03273
03274 gd->work = gdArr[i].work;
03275 gd->bandgap = gdArr[i].bandgap;
03276 gd->therm_eff = gdArr[i].therm_eff;
03277 break;
03278 case STOGNIENKO95:
03279 case BRUGGEMAN35:
03280
03281 gd->work = min(gd->work,gdArr[i].work);
03282 gd->bandgap = min(gd->bandgap,gdArr[i].bandgap);
03283 gd->therm_eff += frac2[i]*gdArr[i].therm_eff;
03284 break;
03285 default:
03286 fprintf( ioQQQ, " Insanity in mie_read_mix\n" );
03287 ShowMe();
03288 cdEXIT(EXIT_FAILURE);
03289 }
03290 gd->matType = gdArr[i].matType;
03291 gd->subl_temp = min(gd->subl_temp,gdArr[i].subl_temp);
03292 }
03293 }
03294
03295 if( gd->rho <= 0. )
03296 {
03297 fprintf( ioQQQ, " Illegal value for the density: %.3e\n", gd->rho );
03298 cdEXIT(EXIT_FAILURE);
03299 }
03300 if( gd->mol_weight <= 0. )
03301 {
03302 fprintf( ioQQQ, " Illegal value for the molecular weight: %.3e\n", gd->mol_weight );
03303 cdEXIT(EXIT_FAILURE);
03304 }
03305 if( maxIndex <= 0. )
03306 {
03307 fprintf( ioQQQ, " No atoms were found in the grain molecule\n" );
03308 cdEXIT(EXIT_FAILURE);
03309 }
03310
03311
03312 ASSERT( wavLo > 0. && wavHi < DBL_MAX && wavLo < wavHi );
03313 ASSERT( gd->abun > 0. && gd->abun < DBL_MAX );
03314 ASSERT( gd->work > 0. && gd->work < DBL_MAX );
03315 ASSERT( gd->bandgap >= 0. && gd->bandgap < gd->work );
03316 ASSERT( gd->therm_eff > 0. && gd->therm_eff <= 1. );
03317 ASSERT( gd->subl_temp > 0. && gd->subl_temp < DBL_MAX );
03318 ASSERT( minIndex > 0. && minIndex < DBL_MAX );
03319
03320
03321 wavLo *= 1. + 10.*DBL_EPSILON;
03322 wavHi *= 1. - 10.*DBL_EPSILON;
03323
03324
03325 nAtoms = 0.;
03326 for( nelem=0; nelem < LIMELM; nelem++ )
03327 {
03328 gd->elmAbun[nelem] /= minIndex;
03329 nAtoms += gd->elmAbun[nelem];
03330 }
03331 ASSERT( nAtoms > 0. );
03332 gd->abun *= minIndex;
03333 gd->mol_weight /= minIndex;
03334
03335 gd->atom_weight = gd->mol_weight/nAtoms;
03336
03337 mie_write_form(gd->elmAbun,chWord);
03338 fprintf( ioQQQ, "\n The chemical formula of the new grain molecule is: %s\n", chWord );
03339 fprintf( ioQQQ, " The abundance wrt H at maximum depletion of this molecule is: %.3e\n",
03340 gd->abun );
03341 fprintf( ioQQQ, " The abundance wrt H at standard depletion of this molecule is: %.3e\n",
03342 gd->abun*gd->depl );
03343
03344
03345 for( nelem=0; nelem < LIMELM; nelem++ )
03346 {
03347 gd->elmAbun[nelem] *= gd->abun*gd->depl;
03348 }
03349
03350 vector<double> delta(sumAxes);
03351 vector<double> frdelta(sumAxes);
03352 vector< complex<double> > eps(sumAxes);
03353
03354 l = 0;
03355 for( i=0; i < nMaterial; i++ )
03356 {
03357 for( k=0; k < gdArr[i].nAxes; k++ )
03358 {
03359 frdelta[l] = gdArr[i].wt[k]*frac[i];
03360 delta[l] = ( l == 0 ) ? frdelta[l] : delta[l-1] + frdelta[l];
03361 ++l;
03362 }
03363 }
03364 ASSERT( l == sumAxes && fabs(delta[l-1]-1.) < 10.*DBL_EPSILON );
03365
03366
03367 gd->wavlen[0].resize( gd->ndata[0] );
03368 gd->n[0].resize( gd->ndata[0] );
03369 gd->nr1[0].resize( gd->ndata[0] );
03370
03371 wavStep = log(wavHi/wavLo)/(double)(gd->ndata[0]-1);
03372
03373 switch( EMTtype )
03374 {
03375 case FARAFONOV00:
03376
03377
03378
03379
03380
03381
03382 for( j=0; j < gd->ndata[0]; j++ )
03383 {
03384 double nre,nim;
03385 complex<double> a1,a2,a1c,a2c,a11,a12,a21,a22,ratio;
03386
03387 gd->wavlen[0][j] = wavLo*exp((double)j*wavStep);
03388
03389 init_eps(gd->wavlen[0][j],nMaterial,gdArr,eps);
03390
03391 ratio = eps[0]/eps[1];
03392
03393 a1 = (ratio+2.)/3.;
03394 a2 = (1.-ratio)*delta[0];
03395
03396 for( l=1; l < sumAxes-1; l++ )
03397 {
03398 ratio = eps[l]/eps[l+1];
03399
03400 a1c = a1;
03401 a2c = a2;
03402 a11 = (ratio+2.)/3.;
03403 a12 = (2.-2.*ratio)/(9.*delta[l]);
03404 a21 = (1.-ratio)*delta[l];
03405 a22 = (2.*ratio+1.)/3.;
03406
03407 a1 = a11*a1c + a12*a2c;
03408 a2 = a21*a1c + a22*a2c;
03409 }
03410
03411 a1c = a1;
03412 a2c = a2;
03413 a11 = 1.;
03414 a12 = 1./3.;
03415 a21 = eps[sumAxes-1];
03416 a22 = -2./3.*eps[sumAxes-1];
03417
03418 a1 = a11*a1c + a12*a2c;
03419 a2 = a21*a1c + a22*a2c;
03420
03421 ratio = a2/a1;
03422 dftori(&nre,&nim,ratio.real(),ratio.imag());
03423
03424 gd->n[0][j] = complex<double>(nre,nim);
03425 gd->nr1[0][j] = nre-1.;
03426 }
03427 break;
03428 case STOGNIENKO95:
03429 case BRUGGEMAN35:
03430 for( j=0; j < gd->ndata[0]; j++ )
03431 {
03432 const double EPS_TOLER = 10.*DBL_EPSILON;
03433 double nre,nim;
03434 complex<double> eps0;
03435
03436 gd->wavlen[0][j] = wavLo*exp((double)j*wavStep);
03437
03438 init_eps(gd->wavlen[0][j],nMaterial,gdArr,eps);
03439
03440
03441 if( j == 0 )
03442 {
03443
03444 eps0 = 0.;
03445 for( l=0; l < sumAxes; l++ )
03446 eps0 += frdelta[l]*eps[l];
03447 }
03448 else
03449 {
03450
03451 eps0 = eps_eff;
03452 }
03453
03454 if( EMTtype == STOGNIENKO95 )
03455
03456
03457 eps_eff = cnewton( Stognienko, frdelta, eps, sumAxes, eps0, EPS_TOLER );
03458 else if( EMTtype == BRUGGEMAN35 )
03459
03460
03461 eps_eff = cnewton( Bruggeman, frdelta, eps, sumAxes, eps0, EPS_TOLER );
03462 else
03463 {
03464 fprintf( ioQQQ, " Insanity in mie_read_mix\n" );
03465 ShowMe();
03466 cdEXIT(EXIT_FAILURE);
03467 }
03468
03469 dftori(&nre,&nim,eps_eff.real(),eps_eff.imag());
03470
03471 gd->n[0][j] = complex<double>(nre,nim);
03472 gd->nr1[0][j] = nre-1.;
03473 }
03474 break;
03475 default:
03476 fprintf( ioQQQ, " Insanity in mie_read_mix\n" );
03477 ShowMe();
03478 cdEXIT(EXIT_FAILURE);
03479 }
03480 return;
03481 }
03482
03483
03484
03485 STATIC void init_eps(double wavlen,
03486 long nMaterial,
03487 const vector<grain_data>& gdArr,
03488 vector< complex<double> >& eps)
03489 {
03490 long i,
03491 k;
03492
03493 long l = 0;
03494
03495 DEBUG_ENTRY( "init_eps()" );
03496 for( i=0; i < nMaterial; i++ )
03497 {
03498 for( k=0; k < gdArr[i].nAxes; k++ )
03499 {
03500 bool lgErr;
03501 long ind;
03502 double eps1,eps2,frc,nim,nre;
03503
03504 find_arr(wavlen,gdArr[i].wavlen[k],gdArr[i].ndata[k],&ind,&lgErr);
03505 ASSERT( !lgErr );
03506 frc = (wavlen-gdArr[i].wavlen[k][ind])/(gdArr[i].wavlen[k][ind+1]-gdArr[i].wavlen[k][ind]);
03507 ASSERT( frc > 0.-10.*DBL_EPSILON && frc < 1.+10.*DBL_EPSILON );
03508 nre = (1.-frc)*gdArr[i].n[k][ind].real() + frc*gdArr[i].n[k][ind+1].real();
03509 ASSERT( nre > 0. );
03510 nim = (1.-frc)*gdArr[i].n[k][ind].imag() + frc*gdArr[i].n[k][ind+1].imag();
03511 ASSERT( nim >= 0. );
03512 ritodf(nre,nim,&eps1,&eps2);
03513 eps[l++] = complex<double>(eps1,eps2);
03514 }
03515 }
03516 return;
03517 }
03518
03519
03520
03521
03522
03523
03524
03525
03526
03527
03528
03529 STATIC complex<double> cnewton(
03530 void(*fun)(complex<double>,const vector<double>&,const vector< complex<double> >&,
03531 long,complex<double>*,double*,double*),
03532 const vector<double>& frdelta,
03533 const vector< complex<double> >& eps,
03534 long sumAxes,
03535 complex<double> x0,
03536 double tol)
03537 {
03538 int i;
03539
03540 const int LOOP_MAX = 100;
03541 const double TINY = 1.e-12;
03542
03543 DEBUG_ENTRY( "cnewton()" );
03544 for( i=0; i < LOOP_MAX; i++ )
03545 {
03546 complex<double> x1,y;
03547 double dudx,dudy,norm2;
03548
03549 (*fun)(x0,frdelta,eps,sumAxes,&y,&dudx,&dudy);
03550
03551 norm2 = pow2(dudx) + pow2(dudy);
03552
03553 if( norm2 < TINY*norm(y) )
03554 {
03555 fprintf( ioQQQ, " cnewton - zero divide error\n" );
03556 ShowMe();
03557 cdEXIT(EXIT_FAILURE);
03558 }
03559 x1 = x0 - complex<double>( y.real()*dudx-y.imag()*dudy, y.imag()*dudx+y.real()*dudy )/norm2;
03560
03561
03562 if( fabs(x0.real()/x1.real()-1.) + fabs(x0.imag()/x1.imag()-1.) < tol )
03563 {
03564 return x1;
03565 }
03566
03567 x0 = x1;
03568 }
03569
03570 fprintf( ioQQQ, " cnewton did not converge\n" );
03571 ShowMe();
03572 cdEXIT(EXIT_FAILURE);
03573 }
03574
03575
03576
03577
03578
03579 STATIC void Stognienko(complex<double> x,
03580 const vector<double>& frdelta,
03581 const vector< complex<double> >& eps,
03582 long sumAxes,
03583 complex<double> *f,
03584 double *dudx,
03585 double *dudy)
03586 {
03587 long i,
03588 l;
03589
03590 static const double L[4] = { 0., 1./2., 1., 1./3. };
03591 static const double fl[4] = { 5./9., 2./9., 2./9., 1. };
03592
03593 DEBUG_ENTRY( "Stognienko()" );
03594 *f = complex<double>(0.,0.);
03595 *dudx = 0.;
03596 *dudy = 0.;
03597 for( l=0; l < sumAxes; l++ )
03598 {
03599 complex<double> hlp = eps[l] - x;
03600 double h1 = eps[l].imag()*x.real() - eps[l].real()*x.imag();
03601
03602 for( i=0; i < 4; i++ )
03603 {
03604 double f1 = fl[i]*frdelta[l];
03605 double xx = ( i < 3 ) ? sin(PI*frdelta[l]) : cos(PI*frdelta[l]);
03606 complex<double> f2 = f1*xx*xx;
03607 complex<double> one = x + hlp*L[i];
03608 complex<double> two = f2*hlp/one;
03609 double h2 = norm(one);
03610 *f += two;
03611 *dudx -= f2.real()*(eps[l].real()*h2 + h1*2.*one.imag()*(1.-L[i]))/pow2(h2);
03612 *dudy -= f2.real()*(eps[l].imag()*h2 - h1*2.*one.real()*(1.-L[i]))/pow2(h2);
03613 }
03614 }
03615 return;
03616 }
03617
03618
03619
03620
03621 STATIC void Bruggeman(complex<double> x,
03622 const vector<double>& frdelta,
03623 const vector< complex<double> >& eps,
03624 long sumAxes,
03625 complex<double> *f,
03626 double *dudx,
03627 double *dudy)
03628 {
03629 long l;
03630
03631 static const double L = 1./3.;
03632
03633 DEBUG_ENTRY( "Bruggeman()" );
03634 *f = complex<double>(0.,0.);
03635 *dudx = 0.;
03636 *dudy = 0.;
03637 for( l=0; l < sumAxes; l++ )
03638 {
03639 complex<double> hlp = eps[l] - x;
03640 double h1 = eps[l].imag()*x.real() - eps[l].real()*x.imag();
03641 complex<double> f2 = frdelta[l];
03642 complex<double> one = x + hlp*L;
03643 complex<double> two = f2*hlp/one;
03644 double h2 = norm(one);
03645 *f += two;
03646 *dudx -= f2.real()*(eps[l].real()*h2 + h1*2.*one.imag()*(1.-L))/pow2(h2);
03647 *dudy -= f2.real()*(eps[l].imag()*h2 - h1*2.*one.real()*(1.-L))/pow2(h2);
03648 }
03649 return;
03650 }
03651
03652
03653
03654 STATIC void mie_read_szd( const char *chFile,
03655 sd_data *sd)
03656 {
03657 bool lgTryOverride = false;
03658 long int dl,
03659 i;
03660 double mul = 0.,
03661 ref_neg = DBL_MAX,
03662 ref_pos = DBL_MAX,
03663 step_neg = DBL_MAX,
03664 step_pos = DBL_MAX;
03665 char chLine[FILENAME_PATH_LENGTH_2],
03666 chWord[WORDLEN];
03667 FILE *io2;
03668
03669
03670
03671
03672
03673
03674 const double FRAC_CUTOFF = 1.e-4;
03675 const double MUL_CO1 = -log(FRAC_CUTOFF);
03676 const double MUL_CO2 = sqrt(MUL_CO1);
03677 const double MUL_CO3 = pow(MUL_CO1,1./3.);
03678 const double MUL_LND = sqrt(-2.*log(FRAC_CUTOFF));
03679 const double MUL_NRM = MUL_LND;
03680
03681 DEBUG_ENTRY( "mie_read_szd()" );
03682
03683 io2 = open_data( chFile, "r", AS_LOCAL_ONLY );
03684
03685 dl = 0;
03686
03687
03688 mie_next_data(chFile,io2,chLine,&dl);
03689 mie_read_long(chFile,chLine,&sd->magic,true,dl);
03690 if( sd->magic != MAGIC_SZD )
03691 {
03692 fprintf( ioQQQ, " Size distribution file %s has obsolete magic number\n",chFile );
03693 fprintf( ioQQQ, " I found %ld, but expected %ld on line #%ld\n",sd->magic,MAGIC_SZD,dl );
03694 fprintf( ioQQQ, " Please replace this file with an up to date version\n" );
03695 cdEXIT(EXIT_FAILURE);
03696 }
03697
03698
03699 mie_next_data(chFile,io2,chLine,&dl);
03700 mie_read_word(chLine,chWord,WORDLEN,true);
03701
03702 if( nMatch( "SSIZ", chWord ) )
03703 {
03704 sd->sdCase = SD_SINGLE_SIZE;
03705 }
03706 else if( nMatch( "NCAR", chWord ) )
03707 {
03708 sd->sdCase = SD_NR_CARBON;
03709 }
03710 else if( nMatch( "POWE", chWord ) )
03711 {
03712 sd->sdCase = SD_POWERLAW;
03713 }
03714 else if( nMatch( "EXP1", chWord ) )
03715 {
03716 sd->sdCase = SD_EXP_CUTOFF1;
03717 sd->a[ipAlpha] = 1.;
03718 mul = MUL_CO1;
03719 }
03720 else if( nMatch( "EXP2", chWord ) )
03721 {
03722 sd->sdCase = SD_EXP_CUTOFF2;
03723 sd->a[ipAlpha] = 2.;
03724 mul = MUL_CO2;
03725 }
03726 else if( nMatch( "EXP3", chWord ) )
03727 {
03728 sd->sdCase = SD_EXP_CUTOFF3;
03729 sd->a[ipAlpha] = 3.;
03730 mul = MUL_CO3;
03731 }
03732 else if( nMatch( "LOGN", chWord ) )
03733 {
03734 sd->sdCase = SD_LOG_NORMAL;
03735 mul = MUL_LND;
03736 }
03737
03738 else if( nMatch( "NORM", chWord ) )
03739 {
03740 sd->sdCase = SD_LIN_NORMAL;
03741 mul = MUL_NRM;
03742 }
03743 else if( nMatch( "TABL", chWord ) )
03744 {
03745 sd->sdCase = SD_TABLE;
03746 }
03747 else
03748 {
03749 sd->sdCase = SD_ILLEGAL;
03750 }
03751
03752 switch( sd->sdCase )
03753 {
03754 case SD_SINGLE_SIZE:
03755
03756 mie_next_data(chFile,io2,chLine,&dl);
03757 mie_read_double(chFile,chLine,&sd->a[ipSize],true,dl);
03758 if( sd->a[ipSize] < SMALLEST_GRAIN || sd->a[ipSize] > LARGEST_GRAIN )
03759 {
03760 fprintf( ioQQQ, " Illegal value for grain size\n" );
03761 fprintf( ioQQQ, " Grain sizes should be between %.5f and %.0f micron\n",
03762 SMALLEST_GRAIN, LARGEST_GRAIN );
03763 fprintf( ioQQQ, " Line #%ld: %s\n",dl,chLine);
03764 cdEXIT(EXIT_FAILURE);
03765 }
03766 break;
03767 case SD_NR_CARBON:
03768
03769 mie_next_data(chFile,io2,chLine,&dl);
03770 mie_read_long(chFile,chLine,&sd->nCarbon,true,dl);
03771 break;
03772 case SD_POWERLAW:
03773
03774 mie_next_data(chFile,io2,chLine,&dl);
03775 mie_read_double(chFile,chLine,&sd->a[ipBLo],true,dl);
03776 if( sd->a[ipBLo] < SMALLEST_GRAIN || sd->a[ipBLo] > LARGEST_GRAIN )
03777 {
03778 fprintf( ioQQQ, " Illegal value for grain size (lower limit)\n" );
03779 fprintf( ioQQQ, " Grain sizes should be between %.5f and %.0f micron\n",
03780 SMALLEST_GRAIN, LARGEST_GRAIN );
03781 fprintf( ioQQQ, " Line #%ld: %s\n",dl,chLine);
03782 cdEXIT(EXIT_FAILURE);
03783 }
03784
03785
03786 mie_next_data(chFile,io2,chLine,&dl);
03787 mie_read_double(chFile,chLine,&sd->a[ipBHi],true,dl);
03788 if( sd->a[ipBHi] < SMALLEST_GRAIN || sd->a[ipBHi] > LARGEST_GRAIN ||
03789 sd->a[ipBHi] <= sd->a[ipBLo] )
03790 {
03791 fprintf( ioQQQ, " Illegal value for grain size (upper limit)\n" );
03792 fprintf( ioQQQ, " Grain sizes should be between %.5f and %.0f micron\n",
03793 SMALLEST_GRAIN, LARGEST_GRAIN );
03794 fprintf( ioQQQ, " and upper limit should be greater than lower limit\n" );
03795 fprintf( ioQQQ, " Line #%ld: %s\n",dl,chLine);
03796 cdEXIT(EXIT_FAILURE);
03797 }
03798
03799
03800 mie_next_data(chFile,io2,chLine,&dl);
03801 if( sscanf( chLine, "%lf", &sd->a[ipExp] ) != 1 )
03802 {
03803 fprintf( ioQQQ, " Syntax error in %s\n",chFile);
03804 fprintf( ioQQQ, " Line #%ld: %s\n",dl,chLine);
03805 cdEXIT(EXIT_FAILURE);
03806 }
03807
03808 sd->a[ipBeta] = 0.;
03809 sd->a[ipSLo] = 0.;
03810 sd->a[ipSHi] = 0.;
03811
03812 sd->lim[ipBLo] = sd->a[ipBLo];
03813 sd->lim[ipBHi] = sd->a[ipBHi];
03814 break;
03815 case SD_EXP_CUTOFF1:
03816 case SD_EXP_CUTOFF2:
03817 case SD_EXP_CUTOFF3:
03818
03819
03820
03821 mie_next_data(chFile,io2,chLine,&dl);
03822 mie_read_double(chFile,chLine,&sd->a[ipBLo],true,dl);
03823
03824
03825 mie_next_data(chFile,io2,chLine,&dl);
03826 mie_read_double(chFile,chLine,&sd->a[ipBHi],true,dl);
03827
03828
03829 mie_next_data(chFile,io2,chLine,&dl);
03830 if( sscanf( chLine, "%lf", &sd->a[ipExp] ) != 1 )
03831 {
03832 fprintf( ioQQQ, " Syntax error in %s\n",chFile);
03833 fprintf( ioQQQ, " Line #%ld: %s\n",dl,chLine);
03834 cdEXIT(EXIT_FAILURE);
03835 }
03836
03837
03838 mie_next_data(chFile,io2,chLine,&dl);
03839 if( sscanf( chLine, "%lf", &sd->a[ipBeta] ) != 1 )
03840 {
03841 fprintf( ioQQQ, " Syntax error in %s\n",chFile);
03842 fprintf( ioQQQ, " Line #%ld: %s\n",dl,chLine);
03843 cdEXIT(EXIT_FAILURE);
03844 }
03845
03846
03847 mie_next_data(chFile,io2,chLine,&dl);
03848 mie_read_double(chFile,chLine,&sd->a[ipSLo],false,dl);
03849
03850
03851 mie_next_data(chFile,io2,chLine,&dl);
03852 mie_read_double(chFile,chLine,&sd->a[ipSHi],false,dl);
03853
03854 ref_neg = sd->a[ipBLo];
03855 step_neg = -mul*sd->a[ipSLo];
03856 ref_pos = sd->a[ipBHi];
03857 step_pos = mul*sd->a[ipSHi];
03858 lgTryOverride = true;
03859 break;
03860 case SD_LOG_NORMAL:
03861
03862 mie_next_data(chFile,io2,chLine,&dl);
03863 mie_read_double(chFile,chLine,&sd->a[ipGCen],true,dl);
03864
03865
03866 mie_next_data(chFile,io2,chLine,&dl);
03867 mie_read_double(chFile,chLine,&sd->a[ipGSig],true,dl);
03868
03869
03870 ref_neg = ref_pos = sd->a[ipGCen]*exp(3.*pow2(sd->a[ipGSig]));
03871 step_neg = sd->a[ipGCen]*(exp(-mul*sd->a[ipGSig]) - 1.);
03872 step_pos = sd->a[ipGCen]*(exp(mul*sd->a[ipGSig]) - 1.);
03873 lgTryOverride = true;
03874 break;
03875 case SD_LIN_NORMAL:
03876
03877 mie_next_data(chFile,io2,chLine,&dl);
03878 mie_read_double(chFile,chLine,&sd->a[ipGCen],true,dl);
03879
03880
03881 mie_next_data(chFile,io2,chLine,&dl);
03882 mie_read_double(chFile,chLine,&sd->a[ipGSig],true,dl);
03883
03884
03885 ref_neg = ref_pos = (sd->a[ipGCen] + sqrt(pow2(sd->a[ipGCen]) + 12.*pow2(sd->a[ipGSig])))/2.;
03886 step_neg = -mul*sd->a[ipGSig];
03887 step_pos = mul*sd->a[ipGSig];
03888 lgTryOverride = true;
03889 break;
03890 case SD_TABLE:
03891
03892 mie_next_data(chFile,io2,chLine,&dl);
03893 mie_read_double(chFile,chLine,&sd->a[ipBLo],true,dl);
03894 if( sd->a[ipBLo] < SMALLEST_GRAIN || sd->a[ipBLo] > LARGEST_GRAIN )
03895 {
03896 fprintf( ioQQQ, " Illegal value for grain size (lower limit)\n" );
03897 fprintf( ioQQQ, " Grain sizes should be between %.5f and %.0f micron\n",
03898 SMALLEST_GRAIN, LARGEST_GRAIN );
03899 fprintf( ioQQQ, " Line #%ld: %s\n",dl,chLine);
03900 cdEXIT(EXIT_FAILURE);
03901 }
03902
03903
03904 mie_next_data(chFile,io2,chLine,&dl);
03905 mie_read_double(chFile,chLine,&sd->a[ipBHi],true,dl);
03906 if( sd->a[ipBHi] < SMALLEST_GRAIN || sd->a[ipBHi] > LARGEST_GRAIN ||
03907 sd->a[ipBHi] <= sd->a[ipBLo] )
03908 {
03909 fprintf( ioQQQ, " Illegal value for grain size (upper limit)\n" );
03910 fprintf( ioQQQ, " Grain sizes should be between %.5f and %.0f micron\n",
03911 SMALLEST_GRAIN, LARGEST_GRAIN );
03912 fprintf( ioQQQ, " and upper limit should be greater than lower limit\n" );
03913 fprintf( ioQQQ, " Line #%ld: %s\n",dl,chLine);
03914 cdEXIT(EXIT_FAILURE);
03915 }
03916
03917
03918 mie_next_data(chFile,io2,chLine,&dl);
03919 mie_read_long(chFile,chLine,&sd->npts,true,dl);
03920 if( sd->npts < 2 )
03921 {
03922 fprintf( ioQQQ, " Illegal value for no. of points in table\n" );
03923 fprintf( ioQQQ, " Line #%ld: %s\n",dl,chLine);
03924 cdEXIT(EXIT_FAILURE);
03925 }
03926
03927
03928 sd->ln_a.resize(sd->npts);
03929 sd->ln_a4dNda.resize(sd->npts);
03930
03931
03932 for( i=0; i < sd->npts; ++i )
03933 {
03934 double help1, help2;
03935
03936 mie_next_data(chFile,io2,chLine,&dl);
03937
03938 if( sscanf( chLine, "%le %le", &help1, &help2 ) != 2 )
03939 {
03940 fprintf( ioQQQ, " Syntax error in %s\n",chFile);
03941 fprintf( ioQQQ, " Line #%ld: %s\n",dl,chLine);
03942 cdEXIT(EXIT_FAILURE);
03943 }
03944
03945 if( help1 <= 0. || help2 <= 0. )
03946 {
03947 fprintf( ioQQQ, " Reading table failed on line #%ld of %s\n",dl,chFile );
03948 fprintf( ioQQQ, " Illegal data value %.6e or %.6e\n", help1, help2 );
03949 cdEXIT(EXIT_FAILURE);
03950 }
03951
03952 sd->ln_a[i] = log(help1);
03953 sd->ln_a4dNda[i] = log(help2);
03954
03955 if( i > 0 && sd->ln_a[i] <= sd->ln_a[i-1] )
03956 {
03957 fprintf( ioQQQ, " Reading table failed on line #%ld of %s\n",dl,chFile );
03958 fprintf( ioQQQ, " Grain radii should be monotonically increasing\n" );
03959 cdEXIT(EXIT_FAILURE);
03960 }
03961 }
03962
03963 sd->lim[ipBLo] = sd->a[ipBLo];
03964 sd->lim[ipBHi] = sd->a[ipBHi];
03965 break;
03966 case SD_ILLEGAL:
03967 default:
03968 fprintf( ioQQQ, " unimplemented case for grain size distribution in file %s\n" , chFile );
03969 fprintf( ioQQQ, " Line #%ld: value %s\n",dl,chWord);
03970 cdEXIT(EXIT_FAILURE);
03971 }
03972
03973
03974
03975
03976
03977
03978
03979 if( lgTryOverride )
03980 {
03981 double help;
03982
03983 mie_next_data(chFile,io2,chLine,&dl);
03984 mie_read_double(chFile,chLine,&help,false,dl);
03985 sd->lim[ipBLo] = ( help <= 0. ) ? search_limit(ref_neg,step_neg,FRAC_CUTOFF,*sd) : help;
03986
03987 mie_next_data(chFile,io2,chLine,&dl);
03988 mie_read_double(chFile,chLine,&help,false,dl);
03989 sd->lim[ipBHi] = ( help <= 0. ) ? search_limit(ref_pos,step_pos,FRAC_CUTOFF,*sd) : help;
03990
03991 if( sd->lim[ipBLo] < SMALLEST_GRAIN || sd->lim[ipBHi] > LARGEST_GRAIN ||
03992 sd->lim[ipBHi] <= sd->lim[ipBLo] )
03993 {
03994 fprintf( ioQQQ, " Illegal size limits: lower %.5f and/or upper %.5f\n",
03995 sd->lim[ipBLo], sd->lim[ipBHi] );
03996 fprintf( ioQQQ, " Grain sizes should be between %.5f and %.0f micron\n",
03997 SMALLEST_GRAIN, LARGEST_GRAIN );
03998 fprintf( ioQQQ, " and upper limit should be greater than lower limit\n" );
03999 fprintf( ioQQQ, " Please alter the size distribution file\n" );
04000 cdEXIT(EXIT_FAILURE);
04001 }
04002 }
04003
04004 fclose(io2);
04005 return;
04006 }
04007
04008
04009 STATIC void mie_read_long( const char *chFile,
04010 const char chLine[],
04011 long int *data,
04012 bool lgZeroIllegal,
04013 long int dl)
04014 {
04015 DEBUG_ENTRY( "mie_read_long()" );
04016 if( sscanf( chLine, "%ld", data ) != 1 )
04017 {
04018 fprintf( ioQQQ, " Syntax error in %s\n",chFile);
04019 fprintf( ioQQQ, " Line #%ld: %s\n",dl,chLine);
04020 cdEXIT(EXIT_FAILURE);
04021 }
04022 if( *data < 0 || (*data == 0 && lgZeroIllegal) )
04023 {
04024 fprintf( ioQQQ, " Illegal data value in %s\n",chFile);
04025 fprintf( ioQQQ, " Line #%ld: %ld\n",dl,*data);
04026 cdEXIT(EXIT_FAILURE);
04027 }
04028 return;
04029 }
04030
04031
04032 STATIC void mie_read_realnum( const char *chFile,
04033 const char chLine[],
04034 realnum *data,
04035 bool lgZeroIllegal,
04036 long int dl)
04037 {
04038 DEBUG_ENTRY( "mie_read_realnum()" );
04039 double help;
04040 if( sscanf( chLine, "%lf", &help ) != 1 )
04041 {
04042 fprintf( ioQQQ, " Syntax error in %s\n",chFile);
04043 fprintf( ioQQQ, " Line #%ld: %s\n",dl,chLine);
04044 cdEXIT(EXIT_FAILURE);
04045 }
04046 *data = (realnum)help;
04047 if( *data < 0. || (*data == 0. && lgZeroIllegal) )
04048 {
04049 fprintf( ioQQQ, " Illegal data value in %s\n",chFile);
04050 fprintf( ioQQQ, " Line #%ld: %14.6e\n",dl,*data);
04051 cdEXIT(EXIT_FAILURE);
04052 }
04053 return;
04054 }
04055
04056
04057 STATIC void mie_read_double( const char *chFile,
04058 const char chLine[],
04059 double *data,
04060 bool lgZeroIllegal,
04061 long int dl)
04062 {
04063 DEBUG_ENTRY( "mie_read_double()" );
04064 if( sscanf( chLine, "%lf", data ) != 1 )
04065 {
04066 fprintf( ioQQQ, " Syntax error in %s\n",chFile);
04067 fprintf( ioQQQ, " Line #%ld: %s\n",dl,chLine);
04068 cdEXIT(EXIT_FAILURE);
04069 }
04070 if( *data < 0. || (*data == 0. && lgZeroIllegal) )
04071 {
04072 fprintf( ioQQQ, " Illegal data value in %s\n",chFile);
04073 fprintf( ioQQQ, " Line #%ld: %14.6e\n",dl,*data);
04074 cdEXIT(EXIT_FAILURE);
04075 }
04076 return;
04077 }
04078
04079
04080 STATIC void mie_read_form( const char chWord[],
04081 double elmAbun[],
04082 double *no_atoms,
04083 double *mol_weight)
04084 {
04085 long int nelem,
04086 len;
04087 char chElmName[3];
04088 double frac;
04089
04090 DEBUG_ENTRY( "mie_read_form()" );
04091
04092 *no_atoms = 0.;
04093 *mol_weight = 0.;
04094 string strWord( chWord );
04095 for( nelem=0; nelem < LIMELM; nelem++ )
04096 {
04097 frac = 0.;
04098 strcpy(chElmName,elementnames.chElementSym[nelem]);
04099 if( chElmName[1] == ' ' )
04100 chElmName[1] = '\0';
04101 string::size_type ptr = strWord.find( chElmName );
04102 if( ptr != string::npos )
04103 {
04104 len = (long)strlen(chElmName);
04105
04106 if( !islower((unsigned char)chWord[ptr+len]) )
04107 {
04108 if( isdigit((unsigned char)chWord[ptr+len]) )
04109 {
04110 sscanf(chWord+ptr+len,"%lf",&frac);
04111 }
04112 else
04113 {
04114 frac = 1.;
04115 }
04116 }
04117 }
04118 elmAbun[nelem] = frac;
04119
04120 if( nelem != ipHYDROGEN )
04121 *no_atoms += frac;
04122 *mol_weight += frac*dense.AtomicWeight[nelem];
04123 }
04124
04125 if( *no_atoms == 0. )
04126 *no_atoms = 1.;
04127 return;
04128 }
04129
04130
04131 STATIC void mie_write_form( const double elmAbun[],
04132 char chWord[])
04133 {
04134 long int len,
04135 nelem;
04136
04137 DEBUG_ENTRY( "mie_write_form()" );
04138
04139 len=0;
04140 for( nelem=0; nelem < LIMELM; nelem++ )
04141 {
04142 if( elmAbun[nelem] > 0. )
04143 {
04144 char chElmName[3];
04145 long index100 = nint(100.*elmAbun[nelem]);
04146
04147 strcpy(chElmName,elementnames.chElementSym[nelem]);
04148 if( chElmName[1] == ' ' )
04149 chElmName[1] = '\0';
04150
04151 if( index100 == 100 )
04152 sprintf( &chWord[len], "%s", chElmName );
04153 else if( index100%100 == 0 )
04154 sprintf( &chWord[len], "%s%ld", chElmName, index100/100 );
04155 else
04156 {
04157 double xIndex = (double)index100/100.;
04158 sprintf( &chWord[len], "%s%.2f", chElmName, xIndex );
04159 }
04160 len = (long)strlen( chWord );
04161 ASSERT( len < FILENAME_PATH_LENGTH_2-25 );
04162 }
04163 }
04164 return;
04165 }
04166
04167
04168 STATIC void mie_read_word( const char chLine[],
04169 char chWord[],
04170 long n,
04171 bool lgToUpper)
04172 {
04173 long ip = 0,
04174 op = 0;
04175
04176 DEBUG_ENTRY( "mie_read_word()" );
04177
04178
04179 while( chLine[ip] == ' ' || chLine[ip] == '"' )
04180 ip++;
04181
04182 while( op < n-1 && chLine[ip] != ' ' && chLine[ip] != '"' )
04183 if( lgToUpper )
04184 chWord[op++] = toupper(chLine[ip++]);
04185 else
04186 chWord[op++] = chLine[ip++];
04187
04188 chWord[op] = '\0';
04189 return;
04190 }
04191
04192
04193 STATIC void mie_next_data( const char *chFile,
04194 FILE *io,
04195 char chLine[],
04196 long int *dl)
04197 {
04198 char *str;
04199
04200 DEBUG_ENTRY( "mie_next_data()" );
04201
04202
04203
04204
04205
04206 chLine[0] = '#';
04207 while( chLine[0] == '#' )
04208 {
04209 mie_next_line(chFile,io,chLine,dl);
04210 }
04211
04212
04213 str = strstr_s(chLine,"#");
04214 if( str != NULL )
04215 *str = '\0';
04216 return;
04217 }
04218
04219
04220
04221 STATIC void mie_next_line( const char *chFile,
04222 FILE *io,
04223 char chLine[],
04224 long int *dl)
04225 {
04226 DEBUG_ENTRY( "mie_next_line()" );
04227
04228 if( read_whole_line( chLine, FILENAME_PATH_LENGTH_2, io ) == NULL )
04229 {
04230 fprintf( ioQQQ, " Could not read from %s\n",chFile);
04231 if( feof(io) )
04232 fprintf( ioQQQ, " EOF reached\n");
04233 fprintf( ioQQQ, " This grain opacity file does not have the expected format.\n");
04234 cdEXIT(EXIT_FAILURE);
04235 }
04236 (*dl)++;
04237 return;
04238 }
04239
04240
04241
04242
04243
04244
04245
04246
04247
04248
04249
04250
04251
04252
04253 void gauss_init(long int nn,
04254 double xbot,
04255 double xtop,
04256 const vector<double>& x,
04257 const vector<double>& a,
04258 vector<double>& rr,
04259 vector<double>& ww)
04260 {
04261 long int i;
04262 double bma,
04263 bpa;
04264
04265 DEBUG_ENTRY( "gauss_init()" );
04266
04267 bpa = (xtop+xbot)/2.;
04268 bma = (xtop-xbot)/2.;
04269
04270 for( i=0; i < nn; i++ )
04271 {
04272 rr[i] = bpa + bma*x[nn-1-i];
04273 ww[i] = bma*a[i];
04274 }
04275 return;
04276 }
04277
04278
04279
04280
04281 void gauss_legendre(long int nn,
04282 vector<double>& x,
04283 vector<double>& a)
04284 {
04285 long int i,
04286 iter,
04287 j;
04288 double cc,
04289 csa,
04290 d,
04291 dp1,
04292 dpn = 0.,
04293 dq,
04294 fj,
04295 fn,
04296 pn,
04297 pn1 = 0.,
04298 q,
04299 xt = 0.;
04300
04301 const double SAFETY = 5.;
04302
04303
04304 DEBUG_ENTRY( "gauss_legendre()" );
04305
04306 if( nn%2 == 1 )
04307 {
04308 fprintf( ioQQQ, " Illegal number of abcissas\n" );
04309 cdEXIT(EXIT_FAILURE);
04310 }
04311
04312 vector<double> c(nn);
04313
04314 fn = (double)nn;
04315 csa = 0.;
04316 cc = 2.;
04317 for( j=1; j < nn; j++ )
04318 {
04319 fj = (double)j;
04320
04321
04322
04323 c[j] = pow2(fj)/((fj-0.5)*(fj+0.5));
04324 cc *= c[j];
04325 }
04326
04327 for( i=0; i < nn/2; i++ )
04328 {
04329 switch( i )
04330 {
04331 case 0:
04332 xt = 1. - 2.78/(4. + pow2(fn));
04333 break;
04334 case 1:
04335 xt = xt - 4.1*(1. + 0.06*(1. - 8./fn))*(1. - xt);
04336 break;
04337 case 2:
04338 xt = xt - 1.67*(1. + 0.22*(1. - 8./fn))*(x[0] - xt);
04339 break;
04340 default:
04341 xt = 3.*(x[i-1] - x[i-2]) + x[i-3];
04342 }
04343 d = 1.;
04344 for( iter=1; (iter < 20) && (fabs(d) > DBL_EPSILON); iter++ )
04345 {
04346
04347
04348
04349 pn1 = 0.5;
04350 pn = xt;
04351 dp1 = 0.;
04352 dpn = 1.;
04353 for( j=1; j < nn; j++ )
04354 {
04355
04356
04357 q = 2.*xt*pn - c[j]*pn1;
04358 dq = 2.*xt*dpn - c[j]*dp1 + 2.*pn;
04359 pn1 = pn;
04360 pn = q;
04361 dp1 = dpn;
04362 dpn = dq;
04363 }
04364 d = pn/dpn;
04365 xt -= d;
04366 }
04367 x[i] = xt;
04368 x[nn-1-i] = -xt;
04369
04370
04371 a[i] = cc/(dpn*2.*pn1);
04372 a[nn-1-i] = a[i];
04373 csa += a[i];
04374 }
04375
04376
04377
04378 if( fabs(1.-csa) > SAFETY*fn*DBL_EPSILON )
04379 {
04380 fprintf( ioQQQ, " gauss_legendre failed to converge: delta = %.4e\n", fabs(1.-csa) );
04381 cdEXIT(EXIT_FAILURE);
04382 }
04383 return;
04384 }
04385
04386
04387
04388
04389
04390
04391 void find_arr(double x,
04392 const vector<double>& xa,
04393 long int n,
04394 long int *ind,
04395 bool *lgOutOfBounds)
04396 {
04397 long int i1,
04398 i2,
04399 i3,
04400 sgn,
04401 sgn2;
04402
04403 DEBUG_ENTRY( "find_arr()" );
04404
04405
04406 if( n < 2 )
04407 {
04408 fprintf( ioQQQ, " Invalid array\n");
04409 cdEXIT(EXIT_FAILURE);
04410 }
04411
04412 i1 = 0;
04413 i3 = n-1;
04414 sgn = sign3(xa[i3]-xa[i1]);
04415 if( sgn == 0 )
04416 {
04417 fprintf( ioQQQ, " Ill-ordered array\n");
04418 cdEXIT(EXIT_FAILURE);
04419 }
04420
04421 *lgOutOfBounds = x < min(xa[0],xa[n-1]) || x > max(xa[0],xa[n-1]);
04422 if( *lgOutOfBounds )
04423 {
04424 *ind = -1;
04425 return;
04426 }
04427
04428 i2 = (n-1)/2;
04429 while( (i3-i1) > 1 )
04430 {
04431 sgn2 = sign3(x-xa[i2]);
04432 if( sgn2 != 0 )
04433 {
04434 if( sgn == sgn2 )
04435 {
04436 i1 = i2;
04437 }
04438 else
04439 {
04440 i3 = i2;
04441 }
04442 i2 = (i1+i3)/2;
04443 }
04444 else
04445 {
04446 *ind = i2;
04447 return;
04448 }
04449 }
04450 *ind = i1;
04451 return;
04452 }
04453
04454
04455
04456
04457
04458
04459
04460
04461
04462
04463
04464
04465
04466
04467
04468
04469
04470
04471
04472
04473
04474
04475
04476
04477
04478
04479
04480
04481 static int NMXLIM = 16000;
04482
04483 STATIC void sinpar(double nre,
04484 double nim,
04485 double x,
04486 double *qext,
04487 double *qphase,
04488 double *qscat,
04489 double *ctbrqs,
04490 double *qback,
04491 long int *iflag)
04492 {
04493 long int n,
04494 nmx1,
04495 nmx2,
04496 nn,
04497 nsqbk;
04498 double ectb,
04499 eqext,
04500 eqpha,
04501 eqscat,
04502 error=0.,
04503 error1=0.,
04504 rx,
04505 t1,
04506 t2,
04507 t3,
04508 t4,
04509 t5,
04510 tx,
04511 x3,
04512 x5=0.,
04513 xcut,
04514 xrd;
04515 complex<double> cdum1,
04516 cdum2,
04517 ci,
04518 eqb,
04519 nc,
04520 nc2,
04521 nc212,
04522 qbck,
04523 rrf,
04524 rrfx,
04525 sman,
04526 sman1,
04527 smbn,
04528 smbn1,
04529 tc1,
04530 tc2,
04531 wn,
04532 wn1,
04533 wn2;
04534
04535 DEBUG_ENTRY( "sinpar()" );
04536
04537 vector< complex<double> > a(NMXLIM);
04538
04539 *iflag = 0;
04540 ci = complex<double>(0.,1.);
04541 nc = complex<double>(nre,-nim);
04542 nc2 = nc*nc;
04543 rrf = 1./nc;
04544 rx = 1./x;
04545 rrfx = rrf*rx;
04546
04547
04548
04549
04550
04551
04552
04553
04554
04555
04556
04557
04558
04559
04560
04561
04562 xrd = exp(log(x)/3.);
04563
04564
04565 t1 = x + 4.*xrd + 3.;
04566
04567
04568
04569 if( !(x <= 8. || x >= 4200.) )
04570 t1 += 0.05*xrd + 3.;
04571 t1 *= 1.01;
04572
04573
04574
04575
04576
04577
04578
04579
04580
04581 t4 = x*sqrt(nre*nre+nim*nim);
04582 nmx1 = nint(max(t1,t4)) + 15;
04583
04584 if( nmx1 < NMXLIM )
04585 {
04586 nmx2 = nint(t1);
04587
04588
04589
04590
04591
04592 if( nmx2 <= 4 )
04593 {
04594 nmx2 = 4;
04595 nmx1 = nint(max(4.,t4)) + 15;
04596 }
04597
04598
04599 a[nmx1] = 0.;
04600
04601
04602
04603
04604 for( n=0; n < nmx1; n++ )
04605 {
04606 nn = nmx1 - n;
04607 a[nn-1] = (double)(nn+1)*rrfx - 1./((double)(nn+1)*rrfx+a[nn]);
04608 }
04609
04610 t1 = cos(x);
04611 t2 = sin(x);
04612 wn2 = complex<double>(t1,-t2);
04613 wn1 = complex<double>(t2,t1);
04614 wn = rx*wn1 - wn2;
04615 tc1 = a[0]*rrf + rx;
04616 tc2 = a[0]*nc + rx;
04617 sman = (tc1*wn.real() - wn1.real())/(tc1*wn - wn1);
04618 smbn = (tc2*wn.real() - wn1.real())/(tc2*wn - wn1);
04619
04620
04621
04622
04623 xcut = 3.e-04;
04624 if( x < xcut )
04625 {
04626 nc212 = (nc2-1.)/(nc2+2.);
04627 x3 = pow3(x);
04628 x5 = x3*pow2(x);
04629
04630 sman = ci*2.*x3*nc212*(1./3.+x*x*0.2*(nc2-2.)/(nc2+2.)) + 4.*x5*x*nc212*nc212/9.;
04631 smbn = ci*x5*(nc2-1.)/45.;
04632 }
04633
04634 sman1 = sman;
04635 smbn1 = smbn;
04636 t1 = 1.5;
04637 sman *= t1;
04638 smbn *= t1;
04639
04640 *qext = 2.*(sman.real() + smbn.real());
04641 *qphase = 2.*(sman.imag() + smbn.imag());
04642 nsqbk = -1;
04643 qbck = -2.*(sman - smbn);
04644 *qscat = (norm(sman) + norm(smbn))/.75;
04645
04646 *ctbrqs = 0.0;
04647 n = 2;
04648
04649
04650 while( true )
04651 {
04652 t1 = 2.*(double)n - 1.;
04653 t3 = 2.*(double)n + 1.;
04654 wn2 = wn1;
04655 wn1 = wn;
04656 wn = t1*rx*wn1 - wn2;
04657 cdum1 = a[n-1];
04658 cdum2 = n*rx;
04659 tc1 = cdum1*rrf + cdum2;
04660 tc2 = cdum1*nc + cdum2;
04661 sman = (tc1*wn.real() - wn1.real())/(tc1*wn - wn1);
04662 smbn = (tc2*wn.real() - wn1.real())/(tc2*wn - wn1);
04663
04664
04665
04666 if( x < xcut && n == 2 )
04667 {
04668
04669 sman = ci*x5*(nc2-1.)/(15.*(2.*nc2+3.));
04670 smbn = 0.;
04671 }
04672
04673 eqext = t3*(sman.real() + smbn.real());
04674 *qext += eqext;
04675 eqpha = t3*(sman.imag() + smbn.imag());
04676 *qphase += eqpha;
04677 nsqbk = -nsqbk;
04678 eqb = t3*(sman - smbn)*(double)nsqbk;
04679 qbck += eqb;
04680 tx = norm(sman) + norm(smbn);
04681 eqscat = t3*tx;
04682 *qscat += eqscat;
04683 t2 = (double)(n - 1);
04684 t5 = (double)n;
04685 t4 = t1/(t5*t2);
04686 t2 = (t2*(t5 + 1.))/t5;
04687 ectb = t2*(sman1.real()*sman.real()+sman1.imag()*sman.imag() + smbn1.real()*smbn.real() +
04688 smbn1.imag()*smbn.imag()) +
04689 t4*(sman1.real()*smbn1.real()+sman1.imag()*smbn1.imag());
04690 *ctbrqs += ectb;
04691
04692
04693
04694 if( tx < 1.e-14 )
04695 {
04696
04697 eqext = fabs(eqext/ *qext);
04698 eqpha = fabs(eqpha/ *qphase);
04699 eqscat = fabs(eqscat/ *qscat);
04700 ectb = ( n == 2 ) ? 0. : fabs(ectb/ *ctbrqs);
04701 eqb = complex<double>( fabs(eqb.real()/qbck.real()), fabs(eqb.imag()/qbck.imag()) );
04702
04703 error = MAX4(eqext,eqpha,eqscat,ectb);
04704
04705 error1 = max(eqb.real(),eqb.imag());
04706 if( error < 1.e-07 && error1 < 1.e-04 )
04707 break;
04708
04709
04710
04711
04712 if( x < xcut )
04713 break;
04714 }
04715
04716 smbn1 = smbn;
04717 sman1 = sman;
04718 n++;
04719 if( n > nmx2 )
04720 {
04721 *iflag = 1;
04722 break;
04723 }
04724 }
04725
04726 t1 = 2.*pow2(rx);
04727 *qext *= t1;
04728 *qphase *= t1;
04729 *qback = norm(qbck)*pow2(rx);
04730 *qscat *= t1;
04731 *ctbrqs *= 2.*t1;
04732 }
04733 else
04734 {
04735 *iflag = 2;
04736 }
04737 return;
04738 }
04739
04740
04741 STATIC void anomal(double x,
04742 double *qext,
04743 double *qabs,
04744 double *qphase,
04745 double *xistar,
04746 double delta,
04747 double beta)
04748 {
04749
04750
04751
04752
04753
04754
04755
04756 double xi,
04757 xii;
04758 complex<double> cbigk,
04759 ci,
04760 cw;
04761
04762 DEBUG_ENTRY( "anomal()" );
04763
04764
04765 xi = 2.*x*delta;
04766 xii = 2.*x*beta;
04767
04768 *xistar = sqrt(pow2(xi)+pow2(xii));
04769
04770 ci = complex<double>(0.,1.);
04771 cw = -complex<double>(xi,xii)*ci;
04772 bigk(cw,&cbigk);
04773 *qext = 4.*cbigk.real();
04774 *qphase = 4.*cbigk.imag();
04775 cw = 2.*xii;
04776 bigk(cw,&cbigk);
04777 *qabs = 2.*cbigk.real();
04778
04779 return;
04780 }
04781
04782
04783 STATIC void bigk(complex<double> cw,
04784 complex<double> *cbigk)
04785 {
04786
04787
04788
04789
04790 DEBUG_ENTRY( "bigk()" );
04791
04792 if( abs(cw) < 1.e-2 )
04793 {
04794
04795
04796
04797 *cbigk = cw*((1./3.)-cw*((1./8.)-cw*((1./30.)-cw*((1./144.)-cw*((1./840.)-cw*(1./5760.))))));
04798 }
04799 else
04800 {
04801 *cbigk = 0.5 + (exp(-cw)*(1.+cw)-1.)/(cw*cw);
04802 }
04803 return;
04804 }
04805
04806
04807
04808 STATIC void ritodf(double nr,
04809 double ni,
04810 double *eps1,
04811 double *eps2)
04812 {
04813
04814 DEBUG_ENTRY( "ritodf()" );
04815
04816 *eps1 = nr*nr - ni*ni;
04817 *eps2 = 2.*nr*ni;
04818 return;
04819 }
04820
04821
04822
04823 STATIC void dftori( double *nr,
04824 double *ni,
04825 double eps1,
04826 double eps2)
04827 {
04828 double eps;
04829
04830 DEBUG_ENTRY( "dftori()" );
04831
04832 eps = sqrt(eps2*eps2+eps1*eps1);
04833 *nr = sqrt((eps+eps1)/2.);
04834 ASSERT( *nr > 0. );
04835
04836
04837
04838 *ni = eps2/(2.*(*nr));
04839 return;
04840 }