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