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