cloudy trunk
Loading...
Searching...
No Matches
grains_mie.cpp
Go to the documentation of this file.
1/* This file is part of Cloudy and is copyright (C)1978-2013 by Gary J. Ferland and
2 * others. For conditions of distribution and use see copyright notice in license.txt */
3#include "cddefines.h"
4#include "continuum.h"
5#include "thirdparty.h"
6#include "elementnames.h"
7#include "physconst.h"
8#include "dense.h"
9#include "called.h"
10#include "version.h"
11#include "grainvar.h"
12#include "rfield.h"
13#include "atmdat.h"
14#include "grains.h"
15
16/*=======================================================*
17 *
18 * Mie code for spherical grains.
19 *
20 * Calculates <pi*a^2*Q_abs>, <pi*a^2*Q_sct>, and (1-<g>)
21 * for arbitrary grain species and size distributions.
22 *
23 * This code is derived from the program cmieuvx.f
24 *
25 * Written by: P.G. Martin (CITA), based on the code described in
26 * >>refer grain physics Hansen, J. E., Travis, L. D. 1974, Space Sci. Rev., 16, 527
27 *
28 * Adapted for Cloudy by Peter A.M. van Hoof (University of Kentucky,
29 * Canadian Institute for Theoretical Astrophysics,
30 * Queen's University of Belfast,
31 * Royal Observatory of Belgium)
32 *
33 *=======================================================*/
34
35
36/* these are the magic numbers for the .rfi, .szd, .opc, and .mix files
37 * the first digit is file type, the rest is date (YYMMDD) */
38static const long MAGIC_RFI = 1030103L;
39static const long MAGIC_SZD = 2010403L;
40static const long MAGIC_OPC = 3100827L;
41static const long MAGIC_MIX = 4030103L;
42
43/* >>chng 02 may 28, by Ryan, moved struct complex to cddefines.h to make it available to entire code. */
44
45/* these are the absolute smallest and largest grain sizes we will
46 * consider (in micron). the lower limit gives a grain with on the
47 * order of one atom in it, so it is physically motivated. the upper
48 * limit comes from the series expansions used in the mie theory,
49 * they will have increasingly more problems converging for larger
50 * grains, so this limit is numerically motivated */
51static const double SMALLEST_GRAIN = 0.0001*(1.-10.*DBL_EPSILON);
52static const double LARGEST_GRAIN = 10.*(1.+10.*DBL_EPSILON);
53
54/* maximum no. of parameters for grain size distribution */
55static const int NSD = 7;
56
57/* these are the indices into the parameter array a[NSD],
58 * NB NB -- the numbers defined below should range from 0 to NSD-1 */
59static const int ipSize = 0;
60static const int ipBLo = 0;
61static const int ipBHi = 1;
62static const int ipExp = 2;
63static const int ipBeta = 3;
64static const int ipSLo = 4;
65static const int ipSHi = 5;
66static const int ipAlpha = 6;
67static const int ipGCen = 2;
68static const int ipGSig = 3;
69
70/* these are the types of refractive index files we recognize */
74
75/* these are the types of EMT's we recognize */
79
80/* these are all the size distribution cases we support */
85
86class sd_data {
87 void p_clear1()
88 {
89 xx.clear();
90 aa.clear();
91 rr.clear();
92 ww.clear();
93 ln_a.clear();
94 ln_a4dNda.clear();
95 }
96public:
97 double a[NSD];
98 double lim[2];
99 double clim[2];
100 vector<double> xx;
101 vector<double> aa;
102 vector<double> rr;
103 vector<double> ww;
104 double unity;
105 double unity_bin;
106 double cSize;
107 double radius;
108 double area;
109 double vol;
110 vector<double> ln_a;
111 vector<double> ln_a4dNda;
113 long int nCarbon;
114 long int magic;
115 long int cPart;
116 long int nPart;
117 long int nmul;
118 long int nn;
119 long int npts;
121 void clear()
122 {
123 p_clear1();
124 }
126 {
127 p_clear1();
128 }
129};
130
131/* maximum no. of principal axes for crystalline grains */
132static const int NAX = 3;
133static const int NDAT = 4;
134
136 void p_clear0()
137 {
138 nAxes = 0;
139 nOpcCols = 0;
140 }
141 void p_clear1()
142 {
143 for( int j=0; j < NAX; j++ )
144 {
145 wavlen[j].clear();
146 n[j].clear();
147 nr1[j].clear();
148 }
149 opcAnu.clear();
150 for( int j=0; j < NDAT; j++ )
151 opcData[j].clear();
152 }
153public:
154 vector<double> wavlen[NAX];
155 vector< complex<double> > n[NAX];
156 vector<double> nr1[NAX];
157 vector<double> opcAnu;
158 vector<double> opcData[NDAT];
159 double wt[NAX];
160 double abun;
161 double depl;
162 double elmAbun[LIMELM];
163 double mol_weight;
164 double atom_weight;
165 double rho;
166 double norm;
167 double work;
168 double bandgap;
169 double therm_eff;
170 double subl_temp;
171 long int magic;
172 long int cAxis;
173 long int nAxes;
174 long int ndata[NAX];
175 long int nOpcCols;
176 long int nOpcData;
177 long int charge;
178 mat_type matType;
180 void clear()
181 {
182 p_clear1();
183 p_clear0();
184 }
186 {
187 p_clear0();
188 }
190 {
191 p_clear1();
192 }
193};
194
195static const int WORDLEN = 5;
196
197/* maximum size for grain type labels */
198static const int LABELSUB1 = 3;
199static const int LABELSUB2 = 5;
200static const int LABELSIZE = LABELSUB1 + LABELSUB2 + 4;
201
202/* this is the number of data points used to set up table of optical constants for a mixed grain */
203static const long MIX_TABLE_SIZE = 2000L;
204
205STATIC void mie_auxiliary(/*@partial@*/sd_data*,/*@in@*/const grain_data*,/*@in@*/const char*);
206STATIC bool mie_auxiliary2(/*@partial@*/vector<int>&,/*@partial@*/multi_arr<double,2>&,
207 /*@partial@*/multi_arr<double,2>&,/*@partial@*/multi_arr<double,2>&,long,long);
208STATIC void mie_integrate(/*@partial@*/sd_data*,double,double,/*@out@*/double*);
209STATIC void mie_cs_size_distr(double,/*@partial@*/sd_data*,/*@in@*/const grain_data*,
210 void(*)(double,/*@in@*/const sd_data*,/*@in@*/const grain_data*,
211 /*@out@*/double*,/*@out@*/double*,/*@out@*/double*,/*@out@*/int*),
212 /*@out@*/double*,/*@out@*/double*,/*@out@*/double*,/*@out@*/int*);
213STATIC void mie_cs(double,/*@in@*/const sd_data*,/*@in@*/const grain_data*,/*@out@*/double*,/*@out@*/double*,
214 /*@out@*/double*,/*@out@*/int*);
215STATIC void ld01_fun(/*@in@*/void(*)(double,const sd_data*,const grain_data[],double*,double*,double*,int*),
216 /*@in@*/double,/*@in@*/double,double,/*@in@*/const sd_data*,/*@in@*/const grain_data[],
217 /*@out@*/double*,/*@out@*/double*,/*@out@*/double*,/*@out@*/int*);
218inline void car1_fun(double,/*@in@*/const sd_data*,/*@in@*/const grain_data[],/*@out@*/double*,/*@out@*/double*,
219 /*@out@*/double*,/*@out@*/int*);
220STATIC void pah1_fun(double,/*@in@*/const sd_data*,/*@in@*/const grain_data*,/*@out@*/double*,/*@out@*/double*,
221 /*@out@*/double*,/*@out@*/int*);
222inline void car2_fun(double,/*@in@*/const sd_data*,/*@in@*/const grain_data[],/*@out@*/double*,/*@out@*/double*,
223 /*@out@*/double*,/*@out@*/int*);
224STATIC void pah2_fun(double,/*@in@*/const sd_data*,/*@in@*/const grain_data*,/*@out@*/double*,/*@out@*/double*,
225 /*@out@*/double*,/*@out@*/int*);
226inline void car3_fun(double,/*@in@*/const sd_data*,/*@in@*/const grain_data[],/*@out@*/double*,/*@out@*/double*,
227 /*@out@*/double*,/*@out@*/int*);
228STATIC void pah3_fun(double,/*@in@*/const sd_data*,/*@in@*/const grain_data*,/*@out@*/double*,/*@out@*/double*,
229 /*@out@*/double*,/*@out@*/int*);
230inline double Drude(double,double,double,double);
231STATIC void tbl_fun(double,/*@in@*/const sd_data*,/*@in@*/const grain_data*,/*@out@*/double*,/*@out@*/double*,
232 /*@out@*/double*,/*@out@*/int*);
233STATIC double size_distr(double,/*@in@*/const sd_data*);
234STATIC double search_limit(double,double,double,sd_data);
235STATIC void mie_calc_ial(/*@in@*/const grain_data*,long,/*@out@*/vector<double>&,/*@in@*/const char*,/*@in@*/bool*);
236STATIC void mie_repair(/*@in@*/const char*,long,int,int,/*@in@*/const double[],double[],/*@in@*/vector<int>&,
237 bool,/*@in@*/bool*);
238STATIC double mie_find_slope(/*@in@*/const double[],/*@in@*/const double[],/*@in@*/const vector<int>&,
239 long,long,int,bool,/*@in@*/bool*);
240STATIC void mie_read_rfi(/*@in@*/const char*,/*@out@*/grain_data*);
241STATIC void mie_read_mix(/*@in@*/const char*,/*@out@*/grain_data*);
242STATIC void init_eps(double,long,/*@in@*/const vector<grain_data>&,/*@out@*/vector< complex<double> >&);
243STATIC complex<double> cnewton(
244 void(*)(complex<double>,const vector<double>&,const vector< complex<double> >&,
245 long,complex<double>*,double*,double*),
246 const vector<double>&,const vector< complex<double> >&,long,complex<double>,double);
247STATIC void Stognienko(complex<double>,const vector<double>&,const vector< complex<double> >&,
248 long,complex<double>*,double*,double*);
249STATIC void Bruggeman(complex<double>,const vector<double>&,const vector< complex<double> >&,
250 long,complex<double>*,double*,double*);
251STATIC void mie_read_szd(/*@in@*/const char*,/*@out@*/sd_data*);
252STATIC void mie_read_long(/*@in@*/const char*,/*@in@*/const char[],/*@out@*/long int*,bool,long int);
253STATIC void mie_read_realnum(/*@in@*/const char*,/*@in@*/const char[],/*@out@*/realnum*,bool,long int);
254STATIC void mie_read_double(/*@in@*/const char*,/*@in@*/const char[],/*@out@*/double*,bool,long int);
255STATIC void mie_read_form(/*@in@*/const char*,/*@out@*/double[],/*@out@*/double*,/*@out@*/double*);
256STATIC void mie_write_form(/*@in@*/const double[],/*@out@*/char[]);
257STATIC void mie_read_word(/*@in@*/const char[],/*@out@*/char[],long,bool);
258STATIC void mie_next_data(/*@in@*/const char*,/*@in@*/FILE*,/*@out@*/char*,/*@in@*/long int*);
259STATIC void mie_next_line(/*@in@*/const char*,/*@in@*/FILE*,/*@out@*/char*,/*@in@*/long int*);
260
261/*=======================================================*/
262/* the following five routines are the core of the Mie code supplied by Peter Martin */
263
264STATIC void sinpar(double,double,double,/*@out@*/double*,/*@out@*/double*,/*@out@*/double*,
265 /*@out@*/double*,/*@out@*/double*,/*@out@*/long*);
266STATIC void anomal(double,/*@out@*/double*,/*@out@*/double*,/*@out@*/double*,/*@out@*/double*,double,double);
267STATIC void bigk(complex<double>,/*@out@*/complex<double>*);
268STATIC void ritodf(double,double,/*@out@*/double*,/*@out@*/double*);
269STATIC void dftori(/*@out@*/double*,/*@out@*/double*,double,double);
270
271
272void mie_write_opc(/*@in@*/ const char *rfi_file,
273 /*@in@*/ const char *szd_file,
274 long int nbin)
275{
276 int Error = 0;
277 bool lgErr,
278 lgErrorOccurred,
279 lgWarning;
280 long int i,
281 nelem,
282 p;
283 double cosb,
284 cs_abs,
285 cs_sct,
286 volfrac,
287 volnorm,
288 wavlen;
289 char chGrainLabel[LABELSIZE+1],
290 ext[3],
292 chFile2[FILENAME_PATH_LENGTH_2],
293 hlp1[LABELSUB1+2],
294 hlp2[LABELSUB2+2],
295 *str,
296 chString[FILENAME_PATH_LENGTH_2];
297 time_t timer;
298 FILE *fdes;
299
300
301 /* no. of logarithmic intervals in table printout of size distribution function */
302 const long NPTS_TABLE = 100L;
303
304 DEBUG_ENTRY( "mie_write_opc()" );
305
306 sd_data sd;
307
308 mie_read_szd( szd_file , &sd );
309
310 sd.nPart = ( sd.sdCase == SD_SINGLE_SIZE || sd.sdCase == SD_NR_CARBON ) ? 1 : nbin;
311 if( sd.nPart <= 0 || sd.nPart >= 100 )
312 {
313 fprintf( ioQQQ, " Illegal number of size distribution bins: %ld\n", sd.nPart );
314 fprintf( ioQQQ, " The number should be between 1 and 99.\n" );
316 }
317 sd.lgLogScale = true;
318
319 vector<grain_data> gdArr(2);
320 grain_data& gd = gdArr[0];
321 grain_data& gd2 = gdArr[1];
322
323 string rfi_string ( rfi_file );
324 if( rfi_string.find( ".rfi" ) != string::npos )
325 {
326 mie_read_rfi( rfi_file , &gd );
327 }
328 else if( rfi_string.find( ".mix" ) != string::npos )
329 {
330 mie_read_mix( rfi_file , &gd );
331 }
332 else
333 {
334 fprintf( ioQQQ, " Refractive index file name %s has wrong extention\n", rfi_file );
335 fprintf( ioQQQ, " It should have extention .rfi or .mix.\n" );
337 }
338
339 if( gd.rfiType == OPC_TABLE && sd.nPart > 1 )
340 {
341 fprintf( ioQQQ, " Illegal number of size distribution bins: %ld\n", sd.nPart );
342 fprintf( ioQQQ, " The number should always be 1 for OPC_TABLE files.\n" );
344 }
345 if( gd.rho <= 0. )
346 {
347 fprintf( ioQQQ, " Illegal value for the specific density: %.4e\n", gd.rho );
349 }
350 if( gd.mol_weight <= 0. )
351 {
352 fprintf( ioQQQ, " Illegal value for the molecular weight: %.4e\n", gd.mol_weight );
354 }
355
356 lgWarning = false;
357
358 /* generate output file name from input file names */
359 strcpy(chFile,rfi_file);
360 str = strstr_s(chFile,".");
361
362 if( str != NULL )
363 *str = '\0';
364
365 strcpy(chFile2,szd_file);
366 str = strstr_s(chFile2,".");
367
368 if( str != NULL )
369 *str = '\0';
370
371 if( sd.sdCase != SD_SINGLE_SIZE && sd.sdCase != SD_NR_CARBON )
372 {
373 sprintf(ext,"%02ld",nbin);
374 strcat(strcat(strcat(strcat(strcat(chFile,"_"),chFile2),"_"),ext),".opc");
375 }
376 else
377 {
378 strcat(strcat(strcat(chFile,"_"),chFile2),".opc");
379 }
380
381 mie_auxiliary(&sd,&gd,"init");
382
383 /* number of protons in plasma per average grain volume */
384 gd.norm = sd.vol*gd.rho/(ATOMIC_MASS_UNIT*gd.mol_weight*gd.abun*gd.depl);
385 volnorm = sd.vol;
386 volfrac = 1.;
387
388 multi_arr<double,2> acs_abs( sd.nPart, rfield.nupper );
389 multi_arr<double,2> acs_sct( acs_abs.clone() );
390 multi_arr<double,2> a1g( acs_abs.clone() );
391 vector<double> inv_att_len( rfield.nupper );
392
393 fprintf( ioQQQ, "\n Starting mie_write_opc, output will be written to %s\n\n", chFile );
394
395 fdes = open_data( chFile, "w", AS_LOCAL_ONLY );
396 lgErr = false;
397
398 (void)time(&timer);
399 lgErr = lgErr || ( fprintf(fdes,"# this file was created by Cloudy %s (%s) on %s",
400 t_version::Inst().chVersion,t_version::Inst().chDate,ctime(&timer)) < 0 );
401 lgErr = lgErr || ( fprintf(fdes,"# ===========================================\n#\n") < 0 );
402 lgErr = lgErr || ( fprintf(fdes,"%12ld # magic number opacity file\n",MAGIC_OPC) < 0 );
403 lgErr = lgErr || ( fprintf(fdes,"%12ld # magic number rfi/mix file\n",gd.magic) < 0 );
404 lgErr = lgErr || ( fprintf(fdes,"%12ld # magic number szd file\n",sd.magic) < 0 );
405
406 /* generate grain label for Cloudy output
407 * adjust LABELSIZE in mie.h when the format defined below is changed ! */
408 strncpy(hlp1,chFile,(size_t)(LABELSUB1+1));
409 hlp1[LABELSUB1+1] = '\0';
410 str = strstr_s(hlp1,"-");
411
412 if( str != NULL )
413 *str = '\0';
414
415 strncpy(hlp2,chFile2,(size_t)(LABELSUB2+1));
416 hlp2[LABELSUB2+1] = '\0';
417 str = strstr_s(hlp2,"-");
418
419 if( str != NULL )
420 *str = '\0';
421
422 strcpy(chGrainLabel," ");
423 if( sd.nPart > 1 )
424 {
425 hlp1[LABELSUB1] = '\0';
426 hlp2[LABELSUB2] = '\0';
427 strcat(strcat(strcat(strcat(chGrainLabel,hlp1),"-"),hlp2),"xx");
428 lgErr = lgErr || ( fprintf(fdes,"%-12.12s # grain type label, xx will be replaced by bin no.\n",
429 chGrainLabel) < 0 );
430 }
431 else
432 {
433 strcat(strcat(strcat(chGrainLabel,hlp1),"-"),hlp2);
434 lgErr = lgErr || ( fprintf(fdes,"%-12.12s # grain type label\n", chGrainLabel) < 0 );
435 }
436
437 lgErr = lgErr || ( fprintf(fdes,"%.6e # specific weight (g/cm^3)\n",gd.rho) < 0 );
438 lgErr = lgErr || ( fprintf(fdes,"%.6e # molecular weight of grain molecule (amu)\n",gd.mol_weight) < 0 );
439 lgErr = lgErr || ( fprintf(fdes,"%.6e # average molecular weight per atom (amu)\n", gd.atom_weight) < 0 );
440 lgErr = lgErr || ( fprintf(fdes,"%.6e # abundance of grain molecule at max depletion\n",gd.abun) < 0 );
441 lgErr = lgErr || ( fprintf(fdes,"%.6e # depletion efficiency\n",gd.depl) < 0 );
442 lgErr = lgErr || ( fprintf(fdes,"%.6e # average grain radius <a^3>/<a^2>, full size distr (cm)\n",
443 3.*sd.vol/sd.area) < 0 );
444 lgErr = lgErr || ( fprintf(fdes,"%.6e # average grain surface area <4pi*a^2>, full size distr (cm^2)\n",
445 sd.area) < 0 );
446 lgErr = lgErr || ( fprintf(fdes,"%.6e # average grain volume <4/3pi*a^3>, full size distr (cm^3)\n",
447 sd.vol) < 0 );
448 lgErr = lgErr || ( fprintf(fdes,"%.6e # total grain radius Int(a) per H, full size distr (cm/H)\n",
449 sd.radius/gd.norm) < 0 );
450 lgErr = lgErr || ( fprintf(fdes,"%.6e # total grain area Int(4pi*a^2) per H, full size distr (cm^2/H)\n",
451 sd.area/gd.norm) < 0 );
452 lgErr = lgErr || ( fprintf(fdes,"%.6e # total grain vol Int(4/3pi*a^3) per H, full size distr (cm^3/H)\n",
453 sd.vol/gd.norm) < 0 );
454 lgErr = lgErr || ( fprintf(fdes,"%.6e # work function (Ryd)\n",gd.work) < 0 );
455 lgErr = lgErr || ( fprintf(fdes,"%.6e # gap between valence and conduction band (Ryd)\n",gd.bandgap) < 0 );
456 lgErr = lgErr || ( fprintf(fdes,"%.6e # efficiency of thermionic emission\n",gd.therm_eff) < 0 );
457 lgErr = lgErr || ( fprintf(fdes,"%.6e # sublimation temperature (K)\n",gd.subl_temp) < 0 );
458 lgErr = lgErr || ( fprintf(fdes,"%12d # material type, 1=carbonaceous, 2=silicate\n",gd.matType) < 0 );
459 lgErr = lgErr || ( fprintf(fdes,"#\n# abundances of constituent elements rel. to hydrogen\n#\n") < 0 );
460
461 for( nelem=0; nelem < LIMELM; nelem++ )
462 {
463 lgErr = lgErr || ( fprintf(fdes,"%.6e # %s\n",gd.elmAbun[nelem],
464 elementnames.chElementSym[nelem]) < 0 );
465 }
466
467 if( sd.sdCase != SD_SINGLE_SIZE && sd.sdCase != SD_NR_CARBON )
468 {
469 lgErr = lgErr || ( fprintf(fdes,"#\n# entire size distribution, amin=%.5f amax=%.5f micron\n",
470 sd.lim[ipBLo],sd.lim[ipBHi]) < 0 );
471 lgErr = lgErr || ( fprintf(fdes,"#\n%.6e # ratio a_max/a_min in each size bin\n",
472 pow(sd.lim[ipBHi]/sd.lim[ipBLo],1./(double)sd.nPart) ) < 0 );
473 lgErr = lgErr || ( fprintf(fdes,"#\n# size distribution function\n#\n") < 0 );
474 lgErr = lgErr || ( fprintf(fdes,"%12ld # number of table entries\n#\n",NPTS_TABLE+1) < 0 );
475 lgErr = lgErr || ( fprintf(fdes,"# ============================\n") < 0 );
476 lgErr = lgErr || ( fprintf(fdes,"# size (micr) a^4*dN/da (cm^3/H)\n#\n") < 0 );
477 for( i=0; i <= NPTS_TABLE; i++ )
478 {
479 double radius, a4dNda;
480 radius = sd.lim[ipBLo]*exp((double)i/(double)NPTS_TABLE*log(sd.lim[ipBHi]/sd.lim[ipBLo]));
481 radius = max(min(radius,sd.lim[ipBHi]),sd.lim[ipBLo]);
482 a4dNda = POW4(radius)*size_distr(radius,&sd)/gd.norm*1.e-12/sd.unity;
483 lgErr = lgErr || ( fprintf(fdes,"%.6e %.6e\n",radius,a4dNda) < 0 );
484 }
485 }
486 else
487 {
488 lgErr = lgErr || ( fprintf(fdes,"#\n") < 0 );
489 lgErr = lgErr || ( fprintf(fdes,"%.6e # a_max/a_min = 1 for single sized grain\n", 1. ) < 0 );
490 lgErr = lgErr || ( fprintf(fdes,"%12ld # no size distribution table\n",0L) < 0 );
491 }
492
493 union {
494 double x;
495 uint32 i[2];
496 } u;
497
498 lgErr = lgErr || ( fprintf(fdes,"#\n") < 0 );
499 lgErr = lgErr || ( fprintf(fdes,"%s # check 1\n",continuum.mesh_md5sum.c_str()) < 0 );
500 u.x = double(rfield.emm);
501 if( cpu.i().big_endian() )
502 lgErr = lgErr || ( fprintf(fdes,"%23.8x %8.8x # check 2\n",u.i[0],u.i[1]) < 0 );
503 else
504 lgErr = lgErr || ( fprintf(fdes,"%23.8x %8.8x # check 2\n",u.i[1],u.i[0]) < 0 );
505 u.x = double(rfield.egamry);
506 if( cpu.i().big_endian() )
507 lgErr = lgErr || ( fprintf(fdes,"%23.8x %8.8x # check 3\n",u.i[0],u.i[1]) < 0 );
508 else
509 lgErr = lgErr || ( fprintf(fdes,"%23.8x %8.8x # check 3\n",u.i[1],u.i[0]) < 0 );
510 u.x = continuum.ResolutionScaleFactor;
511 if( cpu.i().big_endian() )
512 lgErr = lgErr || ( fprintf(fdes,"%23.8x %8.8x # check 4\n",u.i[0],u.i[1]) < 0 );
513 else
514 lgErr = lgErr || ( fprintf(fdes,"%23.8x %8.8x # check 4\n",u.i[1],u.i[0]) < 0 );
515 lgErr = lgErr || ( fprintf(fdes,"%32ld # rfield.nupper\n",rfield.nupper) < 0 );
516 lgErr = lgErr || ( fprintf(fdes,"%32ld # number of size distr. bins\n#\n",sd.nPart) < 0 );
517
518 if( gd.rfiType == OPC_PAH1 )
519 {
520 gd2.clear();
521 mie_read_rfi("graphite.rfi",&gd2);
522 }
523 else if( gd.rfiType == OPC_PAH2N || gd.rfiType == OPC_PAH2C ||
524 gd.rfiType == OPC_PAH3N || gd.rfiType == OPC_PAH3C )
525 {
526 gd2.clear();
527 mie_read_rfi("gdraine.rfi",&gd2);
528 }
529
530 vector<int> ErrorIndex( rfield.nupper );
531
532 for( p=0; p < sd.nPart; p++ )
533 {
534 sd.cPart = p;
535
536 mie_auxiliary(&sd,&gd,"step");
537
538 if( sd.nPart > 1 )
539 {
540 /* >>chng 01 mar 20, creating mie_integrate introduced a change in the normalization
541 * of sd.radius, sd.area, and sd.vol; they now give average quantities for this bin.
542 * gd.norm converts average quantities to integrated quantities per H assuming the
543 * number of grains for the entire size distribution, hence multiplication by frac is
544 * needed to convert to the number of grains in this particular size bin, PvH */
545 double frac = sd.unity_bin/sd.unity;
546 volfrac = sd.vol*frac/volnorm;
547 fprintf( ioQQQ, " Starting size bin %ld, amin=%.5f amax=%.5f micron\n",
548 p+1,sd.clim[ipBLo],sd.clim[ipBHi] );
549 lgErr = lgErr || ( fprintf(fdes,"# size bin %ld, amin=%.5f amax=%.5f micron\n",
550 p+1,sd.clim[ipBLo],sd.clim[ipBHi]) < 0 );
551 lgErr = lgErr || ( fprintf(fdes,"%.6e # average grain ",3.*sd.vol/sd.area) < 0 );
552 lgErr = lgErr || ( fprintf(fdes,"radius <a^3>/<a^2>, this bin (cm)\n") < 0 );
553 lgErr = lgErr || ( fprintf(fdes,"%.6e # average ",sd.area) < 0 );
554 lgErr = lgErr || ( fprintf(fdes,"grain area <4pi*a^2>, this bin (cm^2)\n") < 0 );
555 lgErr = lgErr || ( fprintf(fdes,"%.6e # average ",sd.vol) < 0 );
556 lgErr = lgErr || ( fprintf(fdes,"grain volume <4/3pi*a^3>, this bin (cm^3)\n") < 0 );
557 lgErr = lgErr || ( fprintf(fdes,"%.6e # total grain ",sd.radius*frac/gd.norm) < 0 );
558 lgErr = lgErr || ( fprintf(fdes,"radius Int(a) per H, this bin (cm/H)\n") < 0 );
559 lgErr = lgErr || ( fprintf(fdes,"%.6e # total grain area ",sd.area*frac/gd.norm) < 0 );
560 lgErr = lgErr || ( fprintf(fdes,"Int(4pi*a^2) per H, this bin (cm^2/H)\n") < 0 );
561 lgErr = lgErr || ( fprintf(fdes,"%.6e # total grain volume ",sd.vol*frac/gd.norm) < 0 );
562 lgErr = lgErr || ( fprintf(fdes,"Int(4/3pi*a^3) per H, this bin (cm^3/H)\n#\n") < 0 );
563 }
564
565 lgErrorOccurred = false;
566
567 /* calculate the opacity data */
568 for( i=0; i < rfield.nupper; i++ )
569 {
570 wavlen = WAVNRYD/rfield.anu[i]*1.e4;
571
572 ErrorIndex[i] = 0;
573 acs_abs[p][i] = 0.;
574 acs_sct[p][i] = 0.;
575 a1g[p][i] = 0.;
576
577 switch( gd.rfiType )
578 {
579 case RFI_TABLE:
580 for( gd.cAxis=0; gd.cAxis < gd.nAxes; gd.cAxis++ )
581 {
582 mie_cs_size_distr(wavlen,&sd,&gd,mie_cs,&cs_abs,&cs_sct,&cosb,&Error);
583 ErrorIndex[i] = max(ErrorIndex[i],Error);
584 acs_abs[p][i] += cs_abs*gd.wt[gd.cAxis];
585 acs_sct[p][i] += cs_sct*gd.wt[gd.cAxis];
586 a1g[p][i] += cs_sct*(1.-cosb)*gd.wt[gd.cAxis];
587 }
588 lgErrorOccurred = mie_auxiliary2(ErrorIndex,acs_abs,acs_sct,a1g,p,i);
589 break;
590 case OPC_TABLE:
591 gd.cAxis = 0;
592 mie_cs_size_distr(wavlen,&sd,&gd,tbl_fun,&cs_abs,&cs_sct,&cosb,&Error);
593 ErrorIndex[i] = min(Error,2);
594 lgErrorOccurred = lgErrorOccurred || ( Error > 0 );
595 acs_abs[p][i] = cs_abs*gd.norm;
596 acs_sct[p][i] = cs_sct*gd.norm;
597 a1g[p][i] = 1.-cosb;
598 break;
599 case OPC_GREY:
600 ErrorIndex[i] = 0;
601 acs_abs[p][i] = 1.3121e-23*volfrac*gd.norm;
602 acs_sct[p][i] = 2.6242e-23*volfrac*gd.norm;
603 a1g[p][i] = 1.;
604 break;
605 case OPC_PAH1:
606 gd.cAxis = 0;
607 for( gd2.cAxis=0; gd2.cAxis < gd2.nAxes; gd2.cAxis++ )
608 {
609 mie_cs_size_distr(wavlen,&sd,&gd,car1_fun,&cs_abs,&cs_sct,&cosb,&Error);
610 ErrorIndex[i] = max(ErrorIndex[i],Error);
611 acs_abs[p][i] += cs_abs*gd2.wt[gd2.cAxis];
612 acs_sct[p][i] += 0.1*cs_abs*gd2.wt[gd2.cAxis];
613 a1g[p][i] += 0.1*cs_abs*1.*gd2.wt[gd2.cAxis];
614 }
615 lgErrorOccurred = mie_auxiliary2(ErrorIndex,acs_abs,acs_sct,a1g,p,i);
616 break;
617 case OPC_PAH2N:
618 case OPC_PAH2C:
619 gd.cAxis = 0;
620 // any non-zero charge will do in the OPC_PAH2C case
621 gd.charge = ( gd.rfiType == OPC_PAH2N ) ? 0 : 1;
622 for( gd2.cAxis=0; gd2.cAxis < gd2.nAxes; gd2.cAxis++ )
623 {
624 mie_cs_size_distr(wavlen,&sd,&gd,car2_fun,&cs_abs,&cs_sct,&cosb,&Error);
625 ErrorIndex[i] = max(ErrorIndex[i],Error);
626 acs_abs[p][i] += cs_abs*gd2.wt[gd2.cAxis];
627 acs_sct[p][i] += 0.1*cs_abs*gd2.wt[gd2.cAxis];
628 a1g[p][i] += 0.1*cs_abs*1.*gd2.wt[gd2.cAxis];
629 }
630 lgErrorOccurred = mie_auxiliary2(ErrorIndex,acs_abs,acs_sct,a1g,p,i);
631 break;
632 case OPC_PAH3N:
633 case OPC_PAH3C:
634 gd.cAxis = 0;
635 // any non-zero charge will do in the OPC_PAH3C case
636 gd.charge = ( gd.rfiType == OPC_PAH3N ) ? 0 : 1;
637 for( gd2.cAxis=0; gd2.cAxis < gd2.nAxes; gd2.cAxis++ )
638 {
639 mie_cs_size_distr(wavlen,&sd,&gd,car3_fun,&cs_abs,&cs_sct,&cosb,&Error);
640 ErrorIndex[i] = max(ErrorIndex[i],Error);
641 acs_abs[p][i] += cs_abs*gd2.wt[gd2.cAxis];
642 acs_sct[p][i] += 0.1*cs_abs*gd2.wt[gd2.cAxis];
643 a1g[p][i] += 0.1*cs_abs*1.*gd2.wt[gd2.cAxis];
644 }
645 lgErrorOccurred = mie_auxiliary2(ErrorIndex,acs_abs,acs_sct,a1g,p,i);
646 break;
647 default:
648 fprintf( ioQQQ, " Insanity in mie_write_opc\n" );
649 ShowMe();
651 }
652 }
653
654 /* extrapolate/interpolate for missing data */
655 if( lgErrorOccurred )
656 {
657 strcpy(chString,"absorption cs");
658 mie_repair(chString,rfield.nupper,2,0,rfield.anu,&acs_abs[p][0],ErrorIndex,false,&lgWarning);
659 strcpy(chString,"scattering cs");
660 mie_repair(chString,rfield.nupper,2,1,rfield.anu,&acs_sct[p][0],ErrorIndex,false,&lgWarning);
661 strcpy(chString,"asymmetry parameter");
662 mie_repair(chString,rfield.nupper,1,1,rfield.anu,&a1g[p][0],ErrorIndex,true,&lgWarning);
663 }
664
665 for( i=0; i < rfield.nupper; i++ )
666 {
667 acs_abs[p][i] /= gd.norm;
668 /* >>chng 02 dec 30, do not multiply with (1-g) and write this factor out
669 * separately; this is useful for calculating extinction properties of grains, PvH */
670 acs_sct[p][i] /= gd.norm;
671 }
672 }
673
674 lgErr = lgErr || ( fprintf(fdes,"#\n") < 0 );
675 lgErr = lgErr || ( fprintf(fdes,"# ===========================================\n") < 0 );
676 lgErr = lgErr || ( fprintf(fdes,"# anu (Ryd) abs_cs_01 (cm^2/H) abs_cs_02.....\n#\n") < 0 );
677
678 for( i=0; i < rfield.nupper; i++ )
679 {
680 lgErr = lgErr || ( fprintf(fdes,"%.6e ",rfield.anu[i]) < 0 );
681 for( p=0; p < sd.nPart; p++ )
682 {
683 lgErr = lgErr || ( fprintf(fdes,"%.6e ",acs_abs[p][i]) < 0 );
684 }
685 lgErr = lgErr || ( fprintf(fdes,"\n") < 0 );
686 }
687
688 lgErr = lgErr || ( fprintf(fdes,"#\n") < 0 );
689 lgErr = lgErr || ( fprintf(fdes,"# ===========================================\n") < 0 );
690 lgErr = lgErr || ( fprintf(fdes,"# anu (Ryd) sct_cs_01 (cm^2/H) sct_cs_02.....\n#\n") < 0 );
691
692 for( i=0; i < rfield.nupper; i++ )
693 {
694 lgErr = lgErr || ( fprintf(fdes,"%.6e ",rfield.anu[i]) < 0 );
695 for( p=0; p < sd.nPart; p++ )
696 {
697 lgErr = lgErr || ( fprintf(fdes,"%.6e ",acs_sct[p][i]) < 0 );
698 }
699 lgErr = lgErr || ( fprintf(fdes,"\n") < 0 );
700 }
701
702 lgErr = lgErr || ( fprintf(fdes,"#\n") < 0 );
703 lgErr = lgErr || ( fprintf(fdes,"# ===========================================\n") < 0 );
704 lgErr = lgErr || ( fprintf(fdes,"# anu (Ryd) (1-g)_bin_01 (1-g)_bin_02.....\n#\n") < 0 );
705
706 for( i=0; i < rfield.nupper; i++ )
707 {
708 lgErr = lgErr || ( fprintf(fdes,"%.6e ",rfield.anu[i]) < 0 );
709 for( p=0; p < sd.nPart; p++ )
710 {
711 // cap of 1-g is needed when g is negative...
712 lgErr = lgErr || ( fprintf(fdes,"%.6e ",min(a1g[p][i],1.)) < 0 );
713 }
714 lgErr = lgErr || ( fprintf(fdes,"\n") < 0 );
715 }
716
717 fprintf( ioQQQ, " Starting calculation of inverse attenuation length\n" );
718 strcpy(chString,"inverse attenuation length");
719 if( gd.rfiType != RFI_TABLE )
720 {
721 /* >>chng 02 sep 18, added graphite case for special files like PAH's, PvH */
722 gd2.clear();
723 ial_type icase = gv.which_ial[gd.matType];
724 switch( icase )
725 {
726 case IAL_CAR:
727 mie_read_rfi("graphite.rfi",&gd2);
728 mie_calc_ial(&gd2,rfield.nupper,inv_att_len,chString,&lgWarning);
729 break;
730 case IAL_SIL:
731 mie_read_rfi("silicate.rfi",&gd2);
732 mie_calc_ial(&gd2,rfield.nupper,inv_att_len,chString,&lgWarning);
733 break;
734 default:
735 fprintf( ioQQQ, " mie_write_opc detected unknown ial type: %d\n" , icase );
737 }
738 }
739 else
740 {
741 mie_calc_ial(&gd,rfield.nupper,inv_att_len,chString,&lgWarning);
742 }
743
744 lgErr = lgErr || ( fprintf(fdes,"#\n") < 0 );
745 lgErr = lgErr || ( fprintf(fdes,"# ===========================================\n") < 0 );
746 lgErr = lgErr || ( fprintf(fdes,"# anu (Ryd) inverse attenuation length (cm^-1)\n#\n") < 0 );
747
748 for( i=0; i < rfield.nupper; i++ )
749 {
750 lgErr = lgErr || ( fprintf(fdes,"%.6e %.6e\n",rfield.anu[i],inv_att_len[i]) < 0 );
751 }
752
753 fclose(fdes);
754
755 if( lgErr )
756 {
757 fprintf( ioQQQ, "\n Error writing file: %s\n", chFile );
758 if( remove(chFile) == 0 )
759 {
760 fprintf( ioQQQ, " The file has been removed\n" );
762 }
763 }
764 else
765 {
766 fprintf( ioQQQ, "\n Opacity file %s written succesfully\n\n", chFile );
767 if( lgWarning )
768 {
769 fprintf( ioQQQ, "\n !!! Warnings were detected !!!\n\n" );
770 }
771 }
772 return;
773}
774
775
776STATIC void mie_auxiliary(/*@partial@*/ sd_data *sd,
777 /*@in@*/ const grain_data *gd,
778 /*@in@*/ const char *auxCase)
779{
780 double amin,
781 amax,
782 delta,
783 oldvol,
784 step;
785
786 /* desired relative accuracy of integration over size distribution */
787 const double TOLER = 1.e-3;
788
789 DEBUG_ENTRY( "mie_auxiliary()" );
790 if( strcmp(auxCase,"init") == 0 )
791 {
792 double mass, radius, CpMolecule;
793 /* this is the initial estimate for the multiplier needed to get the
794 * number of abscissas in the gaussian quadrature, the correct value
795 * will be iterated below */
796 sd->nmul = 1;
797
798 /* calculate average grain surface area and volume over size distribution */
799 switch( sd->sdCase )
800 {
801 case SD_SINGLE_SIZE:
802 sd->radius = sd->a[ipSize]*1.e-4;
803 sd->area = 4.*PI*pow2(sd->a[ipSize])*1.e-8;
804 sd->vol = 4./3.*PI*pow3(sd->a[ipSize])*1.e-12;
805 break;
806 case SD_NR_CARBON:
807 if( gd->elmAbun[ipCARBON] == 0. )
808 {
809 fprintf( ioQQQ, "\n This size distribution can only be combined with"
810 " carbonaceous material, bailing out...\n" );
812 }
813 // calculate number of C atoms per grain molecule
814 CpMolecule = gd->elmAbun[ipCARBON]/(gd->abun*gd->depl);
815 // now calculate the mass of the whole grain in gram
816 mass = (double)sd->nCarbon/CpMolecule*gd->mol_weight*ATOMIC_MASS_UNIT;
817 radius = pow(3.*mass/(PI4*gd->rho),1./3.);
818 sd->a[ipSize] = radius*1.e4;
819 sd->radius = radius;
820 sd->area = 4.*PI*pow2(radius);
821 sd->vol = 4./3.*PI*pow3(radius);
822 break;
823 case SD_POWERLAW:
824 case SD_EXP_CUTOFF1:
825 case SD_EXP_CUTOFF2:
826 case SD_EXP_CUTOFF3:
827 case SD_LOG_NORMAL:
828 case SD_LIN_NORMAL:
829 case SD_TABLE:
830 /* set up Gaussian quadrature for entire size range,
831 * first estimate no. of abscissas needed */
832 amin = sd->lgLogScale ? log(sd->lim[ipBLo]) : sd->lim[ipBLo];
833 amax = sd->lgLogScale ? log(sd->lim[ipBHi]) : sd->lim[ipBHi];
834
835 sd->clim[ipBLo] = sd->lim[ipBLo];
836 sd->clim[ipBHi] = sd->lim[ipBHi];
837
838 /* iterate nmul until the integrated volume has converged sufficiently */
839 oldvol= 0.;
840 do
841 {
842 sd->nmul *= 2;
843 mie_integrate(sd,amin,amax,&sd->unity);
844 delta = fabs(sd->vol-oldvol)/sd->vol;
845 oldvol = sd->vol;
846 } while( sd->nmul <= 1024 && delta > TOLER );
847
848 if( delta > TOLER )
849 {
850 fprintf( ioQQQ, " could not converge integration of size distribution\n" );
852 }
853
854 /* we can safely reduce nmul by a factor of 2 and
855 * still reach a relative accuracy of TOLER */
856 sd->nmul /= 2;
857 mie_integrate(sd,amin,amax,&sd->unity);
858 break;
859 case SD_ILLEGAL:
860 default:
861 fprintf( ioQQQ, " insane case for grain size distribution: %d\n" , sd->sdCase );
862 ShowMe();
864 }
865 }
866 else if( strcmp(auxCase,"step") == 0 )
867 {
868 /* calculate average grain surface area and volume over size bin */
869 switch( sd->sdCase )
870 {
871 case SD_SINGLE_SIZE:
872 case SD_NR_CARBON:
873 break;
874 case SD_POWERLAW:
875 case SD_EXP_CUTOFF1:
876 case SD_EXP_CUTOFF2:
877 case SD_EXP_CUTOFF3:
878 case SD_LOG_NORMAL:
879 case SD_LIN_NORMAL:
880 case SD_TABLE:
881 amin = sd->lgLogScale ? log(sd->lim[ipBLo]) : sd->lim[ipBLo];
882 amax = sd->lgLogScale ? log(sd->lim[ipBHi]) : sd->lim[ipBHi];
883 step = (amax - amin)/(double)sd->nPart;
884 amin = amin + (double)sd->cPart*step;
885 amax = min(amax,amin + step);
886
887 sd->clim[ipBLo] = sd->lgLogScale ? exp(amin) : amin;
888 sd->clim[ipBHi] = sd->lgLogScale ? exp(amax) : amax;
889
890 sd->clim[ipBLo] = max(sd->clim[ipBLo],sd->lim[ipBLo]);
891 sd->clim[ipBHi] = min(sd->clim[ipBHi],sd->lim[ipBHi]);
892
893 mie_integrate(sd,amin,amax,&sd->unity_bin);
894
895 break;
896 case SD_ILLEGAL:
897 default:
898 fprintf( ioQQQ, " insane case for grain size distribution: %d\n" , sd->sdCase );
899 ShowMe();
901 }
902 }
903 else
904 {
905 fprintf( ioQQQ, " mie_auxiliary called with insane argument: %s\n", auxCase );
906 ShowMe();
908 }
909 return;
910}
911
912
913STATIC bool mie_auxiliary2(vector<int>& ErrorIndex,
914 multi_arr<double,2>& acs_abs,
915 multi_arr<double,2>& acs_sct,
917 long p,
918 long i)
919{
920 DEBUG_ENTRY( "mie_auxiliary2()" );
921
922 bool lgErrorOccurred = false;
923 if( ErrorIndex[i] > 0 )
924 {
925 ErrorIndex[i] = min(ErrorIndex[i],2);
926 lgErrorOccurred = true;
927 }
928
929 switch( ErrorIndex[i] )
930 {
931 /*lint -e616 */
932 case 2:
933 acs_abs[p][i] = 0.;
934 acs_sct[p][i] = 0.;
935 /*lint -fallthrough */
936 /* controls is supposed to flow to the next case */
937 case 1:
938 a1g[p][i] = 0.;
939 break;
940 /*lint +e616 */
941 case 0:
942 a1g[p][i] /= acs_sct[p][i];
943 break;
944 default:
945 fprintf( ioQQQ, " Insane value for ErrorIndex: %d\n", ErrorIndex[i] );
946 ShowMe();
948 }
949
950 /* sanity checks */
951 if( ErrorIndex[i] < 2 )
952 ASSERT( acs_abs[p][i] > 0. && acs_sct[p][i] > 0. );
953 if( ErrorIndex[i] < 1 )
954 ASSERT( a1g[p][i] > 0. );
955
956 return lgErrorOccurred;
957}
958
959
960STATIC void mie_integrate(/*@partial@*/ sd_data *sd,
961 double amin,
962 double amax,
963 /*@out@*/ double *normalization)
964{
965 long int j;
966 double unity;
967
968 DEBUG_ENTRY( "mie_integrate()" );
969
970 /* set up Gaussian quadrature for size range,
971 * first estimate no. of abscissas needed */
972 sd->nn = sd->nmul*((long)(2.*log(sd->clim[ipBHi]/sd->clim[ipBLo])) + 1);
973 sd->nn = min(max(sd->nn,2*sd->nmul),4096);
974 sd->xx.resize(sd->nn);
975 sd->aa.resize(sd->nn);
976 sd->rr.resize(sd->nn);
977 sd->ww.resize(sd->nn);
978 gauss_legendre(sd->nn,sd->xx,sd->aa);
979 gauss_init(sd->nn,amin,amax,sd->xx,sd->aa,sd->rr,sd->ww);
980
981 /* now integrate surface area and volume */
982 unity = 0.;
983 sd->radius = 0.;
984 sd->area = 0.;
985 sd->vol = 0.;
986
987 for( j=0; j < sd->nn; j++ )
988 {
989 double weight;
990
991 /* use extra factor of size in weights when we use logarithmic mesh */
992 if( sd->lgLogScale )
993 {
994 sd->rr[j] = exp(sd->rr[j]);
995 sd->ww[j] *= sd->rr[j];
996 }
997 weight = sd->ww[j]*size_distr(sd->rr[j],sd);
998 unity += weight;
999 sd->radius += weight*sd->rr[j];
1000 sd->area += weight*pow2(sd->rr[j]);
1001 sd->vol += weight*pow3(sd->rr[j]);
1002 }
1003 *normalization = unity;
1004 sd->radius *= 1.e-4/unity;
1005 sd->area *= 4.*PI*1.e-8/unity;
1006 sd->vol *= 4./3.*PI*1.e-12/unity;
1007 return;
1008}
1009
1010/* read in the *.opc file with opacities and other relevant information */
1011void mie_read_opc(/*@in@*/const char *chFile,
1012 /*@in@*/const GrainPar& gp)
1013{
1014 int res,
1015 lgDefaultQHeat;
1016 long int dl,
1017 help,
1018 i,
1019 nelem,
1020 j,
1021 magic,
1022 nbin,
1023 nup;
1024 size_t nd,
1025 nd2;
1026 realnum RefAbund[LIMELM],
1027 VolTotal;
1028 double anu;
1029 double RadiusRatio;
1030 char chLine[FILENAME_PATH_LENGTH_2],
1031 md5sum[NMD5+1],
1032 *str;
1033 FILE *io2;
1034
1035 /* if a_max/a_min in a single size bin is less than
1036 * RATIO_MAX quantum heating will be turned on by default */
1037 const double RATIO_MAX = pow(100.,1./3.);
1038
1039 DEBUG_ENTRY( "mie_read_opc()" );
1040
1041 io2 = open_data( chFile, "r", AS_DATA_LOCAL );
1042
1043 /* include the name of the file we are reading in the Cloudy output */
1044 strcpy( chLine, " " );
1045 if( strlen(chFile) <= 40 )
1046 {
1047 strncpy( &chLine[0], chFile, strlen(chFile) );
1048 }
1049 else
1050 {
1051 strncpy( &chLine[0], chFile, 37 );
1052 strncpy( &chLine[37], "...", 3 );
1053 }
1054 if( called.lgTalk )
1055 fprintf( ioQQQ, " * >>>> mie_read_opc reading file -- %40s <<<< *\n", chLine );
1056
1057 /* >>chng 02 jan 30, check if file has already been read before, PvH */
1058 for( size_t p=0; p < gv.ReadRecord.size(); ++p )
1059 {
1060 if( gv.ReadRecord[p] == chFile )
1061 {
1062 fprintf( ioQQQ, " File %s has already been read before, was this intended ?\n", chFile );
1063 break;
1064 }
1065 }
1066 gv.ReadRecord.push_back( chFile );
1067
1068 /* allocate memory for first bin */
1069 gv.bin.push_back( new GrainBin );
1070 nd = gv.bin.size()-1;
1071
1072 dl = 0; /* line counter for input file */
1073
1074 /* first read magic numbers */
1075 mie_next_data(chFile,io2,chLine,&dl);
1076 mie_read_long(chFile,chLine,&magic,true,dl);
1077 if( magic != MAGIC_OPC )
1078 {
1079 fprintf( ioQQQ, " Opacity file %s has obsolete magic number\n",chFile );
1080 fprintf( ioQQQ, " I found magic number %ld, but expected %ld on line #%ld\n",
1081 magic,MAGIC_OPC,dl );
1082 fprintf( ioQQQ, " Please recompile this file\n" );
1084 }
1085
1086 /* the following two magic numbers are for information only */
1087 mie_next_data(chFile,io2,chLine,&dl);
1088 mie_read_long(chFile,chLine,&magic,true,dl);
1089
1090 mie_next_data(chFile,io2,chLine,&dl);
1091 mie_read_long(chFile,chLine,&magic,true,dl);
1092
1093 /* the grain scale factor is set equal to the abundances scale factor
1094 * that might have appeared on the grains command. Later, in grains.c,
1095 * it will be further multiplied by gv.GrainMetal, the scale factor that
1096 * appears on the metals & grains command. That command may, or may not,
1097 * have been parsed yet, so can't do it at this stage. */
1098 gv.bin[nd]->dstfactor = (realnum)gp.dep;
1099
1100 /* grain type label */
1101 mie_next_data(chFile,io2,chLine,&dl);
1102 strncpy(gv.bin[nd]->chDstLab,chLine,(size_t)LABELSIZE);
1103 gv.bin[nd]->chDstLab[LABELSIZE] = '\0';
1104
1105 /* specific weight (g/cm^3) */
1106 mie_next_data(chFile,io2,chLine,&dl);
1107 mie_read_realnum(chFile,chLine,&gv.bin[nd]->dustp[0],true,dl);
1108 /* constant needed in the evaluation of the electron escape length */
1109 gv.bin[nd]->eec = pow((double)gv.bin[nd]->dustp[0],-0.85);
1110
1111 /* molecular weight of grain molecule (amu) */
1112 mie_next_data(chFile,io2,chLine,&dl);
1113 mie_read_realnum(chFile,chLine,&gv.bin[nd]->dustp[1],true,dl);
1114
1115 /* average molecular weight per atom (amu) */
1116 mie_next_data(chFile,io2,chLine,&dl);
1117 mie_read_realnum(chFile,chLine,&gv.bin[nd]->atomWeight,true,dl);
1118
1119 /* abundance of grain molecule for max depletion */
1120 mie_next_data(chFile,io2,chLine,&dl);
1121 mie_read_realnum(chFile,chLine,&gv.bin[nd]->dustp[2],true,dl);
1122
1123 /* depletion efficiency */
1124 mie_next_data(chFile,io2,chLine,&dl);
1125 mie_read_realnum(chFile,chLine,&gv.bin[nd]->dustp[3],true,dl);
1126
1127 /* fraction of the integrated volume contained in this bin */
1128 gv.bin[nd]->dustp[4] = 1.;
1129
1130 /* average grain radius <a^3>/<a^2> for entire size distr (cm) */
1131 mie_next_data(chFile,io2,chLine,&dl);
1132 mie_read_realnum(chFile,chLine,&gv.bin[nd]->AvRadius,true,dl);
1133 gv.bin[nd]->eyc = 1./gv.bin[nd]->AvRadius + 1.e7;
1134
1135 /* average grain area <4pi*a^2> for entire size distr (cm^2) */
1136 mie_next_data(chFile,io2,chLine,&dl);
1137 mie_read_realnum(chFile,chLine,&gv.bin[nd]->AvArea,true,dl);
1138
1139 /* average grain volume <4/3pi*a^3> for entire size distr (cm^3) */
1140 mie_next_data(chFile,io2,chLine,&dl);
1141 mie_read_realnum(chFile,chLine,&gv.bin[nd]->AvVol,true,dl);
1142
1143 /* total grain radius Int(a) per H for entire size distr (cm/H) */
1144 mie_next_data(chFile,io2,chLine,&dl);
1145 mie_read_realnum(chFile,chLine,&gv.bin[nd]->IntRadius,true,dl);
1146
1147 /* total grain area Int(4pi*a^2) per H for entire size distr (cm^2/H) */
1148 mie_next_data(chFile,io2,chLine,&dl);
1149 mie_read_realnum(chFile,chLine,&gv.bin[nd]->IntArea,true,dl);
1150
1151 /* total grain vol Int(4/3pi*a^3) per H for entire size distr (cm^3/H) */
1152 mie_next_data(chFile,io2,chLine,&dl);
1153 mie_read_realnum(chFile,chLine,&gv.bin[nd]->IntVol,true,dl);
1154
1155 /* work function, in Rydberg */
1156 mie_next_data(chFile,io2,chLine,&dl);
1157 mie_read_realnum(chFile,chLine,&gv.bin[nd]->DustWorkFcn,true,dl);
1158
1159 /* bandgap, in Rydberg */
1160 mie_next_data(chFile,io2,chLine,&dl);
1161 mie_read_realnum(chFile,chLine,&gv.bin[nd]->BandGap,false,dl);
1162
1163 /* efficiency of thermionic emissions, between 0 and 1 */
1164 mie_next_data(chFile,io2,chLine,&dl);
1165 mie_read_realnum(chFile,chLine,&gv.bin[nd]->ThermEff,true,dl);
1166
1167 /* sublimation temperature in K */
1168 mie_next_data(chFile,io2,chLine,&dl);
1169 mie_read_realnum(chFile,chLine,&gv.bin[nd]->Tsublimat,true,dl);
1170
1171 /* material type, determines enthalpy function, etc... */
1172 mie_next_data(chFile,io2,chLine,&dl);
1173 mie_read_long(chFile,chLine,&help,true,dl);
1174 gv.bin[nd]->matType = (mat_type)help;
1175
1176 for( nelem=0; nelem < LIMELM; nelem++ )
1177 {
1178 mie_next_data(chFile,io2,chLine,&dl);
1179 mie_read_realnum(chFile,chLine,&RefAbund[nelem],false,dl);
1180
1181 gv.bin[nd]->elmAbund[nelem] = RefAbund[nelem];
1182
1183 /* this coefficient is defined at the end of appendix A.10 of BFM */
1184 gv.bin[nd]->AccomCoef[nelem] = 2.*gv.bin[nd]->atomWeight*dense.AtomicWeight[nelem]/
1185 pow2(gv.bin[nd]->atomWeight+dense.AtomicWeight[nelem]);
1186 }
1187
1188 /* ratio a_max/a_min for grains in a single size bin */
1189 mie_next_data(chFile,io2,chLine,&dl);
1190 mie_read_double(chFile,chLine,&RadiusRatio,true,dl);
1191
1192 gv.bin[nd]->nDustFunc = gp.nDustFunc;
1193 lgDefaultQHeat = ( RadiusRatio < RATIO_MAX && !gp.lgGreyGrain );
1194 gv.bin[nd]->lgQHeat = ( gp.lgForbidQHeating ) ? false : ( gp.lgRequestQHeating || lgDefaultQHeat );
1195 gv.bin[nd]->cnv_H_pGR = gv.bin[nd]->AvVol/gv.bin[nd]->IntVol;
1196 gv.bin[nd]->cnv_GR_pH = 1./gv.bin[nd]->cnv_H_pGR;
1197
1198 /* this is capacity per grain, in Farad per grain */
1199 gv.bin[nd]->Capacity = PI4*ELECTRIC_CONST*gv.bin[nd]->IntRadius/100.*gv.bin[nd]->cnv_H_pGR;
1200
1201 /* skip the table of the size distribution function (if present) */
1202 mie_next_data(chFile,io2,chLine,&dl);
1203 mie_read_long(chFile,chLine,&nup,false,dl);
1204 for( i=0; i < nup; i++ )
1205 mie_next_data(chFile,io2,chLine,&dl);
1206
1207 /* read checksum of continuum_mesh.ini */
1208 mie_next_data(chFile,io2,chLine,&dl);
1209 mie_read_word(chLine,md5sum,sizeof(md5sum),false);
1210
1211 union {
1212 double x;
1213 uint32 i[2];
1214 } u;
1215 double mesh_lo, mesh_hi;
1216
1217 /* read lower limit of frequency mesh in hex form */
1218 mie_next_data(chFile,io2,chLine,&dl);
1219 if( cpu.i().big_endian() )
1220 sscanf( chLine, "%x %x", &u.i[0], &u.i[1] );
1221 else
1222 sscanf( chLine, "%x %x", &u.i[1], &u.i[0] );
1223 mesh_lo = u.x;
1224
1225 /* read upper limit of frequency mesh in hex form */
1226 mie_next_data(chFile,io2,chLine,&dl);
1227 if( cpu.i().big_endian() )
1228 sscanf( chLine, "%x %x", &u.i[0], &u.i[1] );
1229 else
1230 sscanf( chLine, "%x %x", &u.i[1], &u.i[0] );
1231 mesh_hi = u.x;
1232
1233 if( strncmp( md5sum, continuum.mesh_md5sum.c_str(), NMD5 ) != 0 ||
1234 !fp_equal( mesh_lo, double(rfield.emm) ) ||
1235 !fp_equal( mesh_hi, double(rfield.egamry) ) )
1236 {
1237 fprintf( ioQQQ, " Opacity file %s has an incompatible energy grid.\n", chFile );
1238 fprintf( ioQQQ, " Please recompile this file using the COMPILE GRAINS command.\n" );
1240 }
1241
1242 /* read mesh resolution scale factor in hex form */
1243 mie_next_data(chFile,io2,chLine,&dl);
1244 if( cpu.i().big_endian() )
1245 sscanf( chLine, "%x %x", &u.i[0], &u.i[1] );
1246 else
1247 sscanf( chLine, "%x %x", &u.i[1], &u.i[0] );
1248 /* this number is checked later since it may not have been set yet by the input script */
1249 gv.bin[nd]->RSFCheck = u.x;
1250
1251 /* nup is number of frequency bins stored in file, this should match rfield.nupper */
1252 mie_next_data(chFile,io2,chLine,&dl);
1253 mie_read_long(chFile,chLine,&nup,true,dl);
1254
1255 /* no. of size distribution bins */
1256 mie_next_data(chFile,io2,chLine,&dl);
1257 mie_read_long(chFile,chLine,&nbin,true,dl);
1258
1259 /* now update the fields for a resolved size distribution */
1260 if( nbin > 1 )
1261 {
1262 /* remember this number since it will be overwritten below */
1263 VolTotal = gv.bin[nd]->IntVol;
1264
1265 for( i=0; i < nbin; i++ )
1266 {
1267 if( i == 0 )
1268 nd2 = nd;
1269 else
1270 {
1271 /* allocate memory for remaining bins */
1272 gv.bin.push_back( new GrainBin );
1273 nd2 = gv.bin.size()-1;
1274
1275 /* first do a straight copy of all the fields ... */
1276 *gv.bin[nd2] = *gv.bin[nd];
1277 }
1278
1279 /* ... then update anything that needs updating */
1280
1281 /* average grain radius <a^3>/<a^2> for this bin (cm) */
1282 mie_next_data(chFile,io2,chLine,&dl);
1283 mie_read_realnum(chFile,chLine,&gv.bin[nd2]->AvRadius,true,dl);
1284
1285 /* average grain area in this bin (cm^2) */
1286 mie_next_data(chFile,io2,chLine,&dl);
1287 mie_read_realnum(chFile,chLine,&gv.bin[nd2]->AvArea,true,dl);
1288
1289 /* average grain volume in this bin (cm^3) */
1290 mie_next_data(chFile,io2,chLine,&dl);
1291 mie_read_realnum(chFile,chLine,&gv.bin[nd2]->AvVol,true,dl);
1292
1293 /* total grain radius Int(a) per H for this bin (cm/H) */
1294 mie_next_data(chFile,io2,chLine,&dl);
1295 mie_read_realnum(chFile,chLine,&gv.bin[nd2]->IntRadius,true,dl);
1296
1297 /* total grain area Int(4pi*a^2) per H for this bin (cm^2/H) */
1298 mie_next_data(chFile,io2,chLine,&dl);
1299 mie_read_realnum(chFile,chLine,&gv.bin[nd2]->IntArea,true,dl);
1300
1301 /* total grain vol Int(4/3pi*a^3) per H for this bin (cm^3/H) */
1302 mie_next_data(chFile,io2,chLine,&dl);
1303 mie_read_realnum(chFile,chLine,&gv.bin[nd2]->IntVol,true,dl);
1304
1305 gv.bin[nd2]->cnv_H_pGR = gv.bin[nd2]->AvVol/gv.bin[nd2]->IntVol;
1306 gv.bin[nd2]->cnv_GR_pH = 1./gv.bin[nd2]->cnv_H_pGR;
1307
1308 /* this is capacity per grain, in Farad per grain */
1309 gv.bin[nd2]->Capacity =
1310 PI4*ELECTRIC_CONST*gv.bin[nd2]->IntRadius/100.*gv.bin[nd2]->cnv_H_pGR;
1311
1312 /* dustp[4] gives the fraction of the grain abundance that is
1313 * contained in a particular bin. for unresolved distributions it is
1314 * by definition 1, for resolved distributions it is smaller than 1. */
1315 gv.bin[nd2]->dustp[4] = gv.bin[nd2]->IntVol/VolTotal;
1316 for( nelem=0; nelem < LIMELM; nelem++ )
1317 gv.bin[nd2]->elmAbund[nelem] = RefAbund[nelem]*gv.bin[nd2]->dustp[4];
1318 }
1319
1320 /* this must be in a separate loop! */
1321 for( i=0; i < nbin; i++ )
1322 {
1323 nd2 = nd + i;
1324 /* modify grain labels */
1325 str = strstr_s(gv.bin[nd2]->chDstLab,"xx");
1326 if( str != NULL )
1327 sprintf(str,"%02ld",i+1);
1328 }
1329 }
1330
1331 /* allocate memory for arrays */
1332 for( i=0; i < nbin; i++ )
1333 {
1334 nd2 = nd + i;
1335 gv.bin[nd2]->dstab1.resize(nup);
1336 gv.bin[nd2]->pure_sc1.resize(nup);
1337 gv.bin[nd2]->asym.resize(nup);
1338 gv.bin[nd2]->inv_att_len.resize(nup);
1339 }
1340
1341 /* skip the next 5 lines */
1342 for( i=0; i < 5; i++ )
1343 mie_next_line(chFile,io2,chLine,&dl);
1344
1345 /* now read absorption opacities */
1346 for( i=0; i < nup; i++ )
1347 {
1348 /* read in energy scale and then opacities */
1349 if( (res = fscanf(io2,"%le",&anu)) != 1 )
1350 {
1351 fprintf( ioQQQ, " Read failed on %s\n",chFile );
1352 if( res == EOF )
1353 fprintf( ioQQQ, " EOF reached prematurely\n" );
1355 }
1356 for( j=0; j < nbin; j++ )
1357 {
1358 nd2 = nd + j;
1359 if( (res = fscanf(io2,"%le",&gv.bin[nd2]->dstab1[i])) != 1 )
1360 {
1361 fprintf( ioQQQ, " Read failed on %s\n",chFile );
1362 if( res == EOF )
1363 fprintf( ioQQQ, " EOF reached prematurely\n" );
1365 }
1366 ASSERT( gv.bin[nd2]->dstab1[i] > 0. );
1367 }
1368 }
1369
1370 /* skip to end-of-line and then skip next 4 lines */
1371 for( i=0; i < 5; i++ )
1372 mie_next_line(chFile,io2,chLine,&dl);
1373
1374 /* now read scattering opacities */
1375 for( i=0; i < nup; i++ )
1376 {
1377 if( (res = fscanf(io2,"%le",&anu)) != 1 )
1378 {
1379 fprintf( ioQQQ, " Read failed on %s\n",chFile );
1380 if( res == EOF )
1381 fprintf( ioQQQ, " EOF reached prematurely\n" );
1383 }
1384 for( j=0; j < nbin; j++ )
1385 {
1386 nd2 = nd + j;
1387 if( (res = fscanf(io2,"%le",&gv.bin[nd2]->pure_sc1[i])) != 1 )
1388 {
1389 fprintf( ioQQQ, " Read failed on %s\n",chFile );
1390 if( res == EOF )
1391 fprintf( ioQQQ, " EOF reached prematurely\n" );
1393 }
1394 ASSERT( gv.bin[nd2]->pure_sc1[i] > 0. );
1395 }
1396 }
1397
1398 /* skip to end-of-line and then skip next 4 lines */
1399 for( i=0; i < 5; i++ )
1400 mie_next_line(chFile,io2,chLine,&dl);
1401
1402 /* now read asymmetry factor */
1403 for( i=0; i < nup; i++ )
1404 {
1405 if( (res = fscanf(io2,"%le",&anu)) != 1 )
1406 {
1407 fprintf( ioQQQ, " Read failed on %s\n",chFile );
1408 if( res == EOF )
1409 fprintf( ioQQQ, " EOF reached prematurely\n" );
1411 }
1412 for( j=0; j < nbin; j++ )
1413 {
1414 nd2 = nd + j;
1415 if( (res = fscanf(io2,"%le",&gv.bin[nd2]->asym[i])) != 1 )
1416 {
1417 fprintf( ioQQQ, " Read failed on %s\n",chFile );
1418 if( res == EOF )
1419 fprintf( ioQQQ, " EOF reached prematurely\n" );
1421 }
1422 ASSERT( gv.bin[nd2]->asym[i] > 0. );
1423 // just in case we read an old opacity file...
1424 gv.bin[nd2]->asym[i] = min(gv.bin[nd2]->asym[i],1.);
1425 }
1426 }
1427
1428 /* skip to end-of-line and then skip next 4 lines */
1429 for( i=0; i < 5; i++ )
1430 mie_next_line(chFile,io2,chLine,&dl);
1431
1432 /* now read inverse attenuation length */
1433 for( i=0; i < nup; i++ )
1434 {
1435 double help;
1436 if( (res = fscanf(io2,"%le %le",&anu,&help)) != 2 )
1437 {
1438 fprintf( ioQQQ, " Read failed on %s\n",chFile );
1439 if( res == EOF )
1440 fprintf( ioQQQ, " EOF reached prematurely\n" );
1442 }
1443 gv.bin[nd]->inv_att_len[i] = (realnum)help;
1444 ASSERT( gv.bin[nd]->inv_att_len[i] > 0. );
1445
1446 for( j=1; j < nbin; j++ )
1447 {
1448 nd2 = nd + j;
1449 gv.bin[nd2]->inv_att_len[i] = gv.bin[nd]->inv_att_len[i];
1450 }
1451 }
1452
1453 fclose(io2);
1454 return;
1455}
1456
1457
1458/* calculate average absorption, scattering cross section (i.e. pi a^2 Q) and
1459 * average asymmetry parameter g for an arbitrary grain size distribution */
1460STATIC void mie_cs_size_distr(double wavlen, /* micron */
1461 /*@partial@*/ sd_data *sd,
1462 /*@in@*/ const grain_data *gd,
1463 void(*cs_fun)(double,/*@in@*/const sd_data*,
1464 /*@in@*/const grain_data*,
1465 /*@out@*/double*,/*@out@*/double*,
1466 /*@out@*/double*,/*@out@*/int*),
1467 /*@out@*/ double *cs_abs, /* cm^2, average */
1468 /*@out@*/ double *cs_sct, /* cm^2, average */
1469 /*@out@*/ double *cosb,
1470 /*@out@*/ int *error)
1471{
1472 bool lgADLused;
1473 long int i;
1474 double absval,
1475 g,
1476 sctval,
1477 weight;
1478
1479 DEBUG_ENTRY( "mie_cs_size_distr()" );
1480
1481 /* sanity checks */
1482 ASSERT( wavlen > 0. );
1483 ASSERT( gd->cAxis >= 0 && gd->cAxis < gd->nAxes && gd->cAxis < NAX );
1484
1485 switch( sd->sdCase )
1486 {
1487 case SD_SINGLE_SIZE:
1488 case SD_NR_CARBON:
1489 /* do single sized grain */
1490 ASSERT( sd->a[ipSize] > 0. );
1491 sd->cSize = sd->a[ipSize];
1492 (*cs_fun)(wavlen,sd,gd,cs_abs,cs_sct,cosb,error);
1493 break;
1494 case SD_POWERLAW:
1495 /* simple powerlaw distribution */
1496 case SD_EXP_CUTOFF1:
1497 case SD_EXP_CUTOFF2:
1498 case SD_EXP_CUTOFF3:
1499 /* powerlaw distribution with exponential cutoff */
1500 case SD_LOG_NORMAL:
1501 /* gaussian distribution in ln(a) */
1502 case SD_LIN_NORMAL:
1503 /* gaussian distribution in a */
1504 case SD_TABLE:
1505 /* user supplied table of a^4*dN/da */
1506 ASSERT( sd->lim[ipBLo] > 0. && sd->lim[ipBHi] > 0. && sd->lim[ipBHi] > sd->lim[ipBLo] );
1507 lgADLused = false;
1508 *cs_abs = 0.;
1509 *cs_sct = 0.;
1510 *cosb = 0.;
1511 for( i=0; i < sd->nn; i++ )
1512 {
1513 sd->cSize = sd->rr[i];
1514 (*cs_fun)(wavlen,sd,gd,&absval,&sctval,&g,error);
1515 if( *error >= 2 )
1516 {
1517 /* mie_cs failed to converge -> integration is invalid */
1518 *cs_abs = -1.;
1519 *cs_sct = -1.;
1520 *cosb = -2.;
1521 return;
1522 }
1523 else if( *error == 1 )
1524 {
1525 /* anomalous diffraction limit used -> g is not valid */
1526 lgADLused = true;
1527 }
1528 weight = sd->ww[i]*size_distr(sd->rr[i],sd);
1529 *cs_abs += weight*absval;
1530 *cs_sct += weight*sctval;
1531 *cosb += weight*sctval*g;
1532 }
1533 if( lgADLused )
1534 {
1535 *error = 1;
1536 *cosb = -2.;
1537 }
1538 else
1539 {
1540 *error = 0;
1541 *cosb /= *cs_sct;
1542 }
1543 *cs_abs /= sd->unity;
1544 *cs_sct /= sd->unity;
1545 break;
1546 case SD_ILLEGAL:
1547 default:
1548 fprintf( ioQQQ, " insane case for grain size distribution: %d\n" , sd->sdCase );
1549 ShowMe();
1551 }
1552 /* sanity checks */
1553 if( *error < 2 )
1554 ASSERT( *cs_abs > 0. && *cs_sct > 0. );
1555 if( *error < 1 )
1556 ASSERT( fabs(*cosb) <= 1.+10.*DBL_EPSILON );
1557 return;
1558}
1559
1560
1561/* calculate absorption, scattering cross section (i.e. pi a^2 Q) and
1562 * asymmetry parameter g (=cosb) for a single sized grain defined by gd */
1563STATIC void mie_cs(double wavlen, /* micron */
1564 /*@in@*/ const sd_data *sd,
1565 /*@in@*/ const grain_data *gd,
1566 /*@out@*/ double *cs_abs, /* cm^2 */
1567 /*@out@*/ double *cs_sct, /* cm^2 */
1568 /*@out@*/ double *cosb,
1569 /*@out@*/ int *error)
1570{
1571 bool lgOutOfBounds;
1572 long int iflag,
1573 ind;
1574 double area,
1575 aqabs,
1576 aqext,
1577 aqphas,
1578 beta,
1579 ctbrqs = -DBL_MAX,
1580 delta,
1581 frac,
1582 nim,
1583 nre,
1584 nr1,
1585 qback,
1586 qext = -DBL_MAX,
1587 qphase,
1588 qscatt = -DBL_MAX,
1589 x,
1590 xistar;
1591
1592 DEBUG_ENTRY( "mie_cs()" );
1593
1594 /* sanity checks, should already have been checked further upstream */
1595 ASSERT( wavlen > 0. );
1596 ASSERT( sd->cSize > 0. );
1597 ASSERT( gd->ndata[gd->cAxis] > 1 && (long)gd->wavlen[gd->cAxis].size() == gd->ndata[gd->cAxis] );
1598
1599 /* first interpolate optical constants */
1600 /* >>chng 02 oct 22, moved calculation of optical constants from mie_cs_size_distr to mie_cs,
1601 * this increases overhead, but makes the code in mie_cs_size_distr more transparent, PvH */
1602 find_arr(wavlen,gd->wavlen[gd->cAxis],gd->ndata[gd->cAxis],&ind,&lgOutOfBounds);
1603
1604 if( lgOutOfBounds )
1605 {
1606 *error = 3;
1607 *cs_abs = -1.;
1608 *cs_sct = -1.;
1609 *cosb = -2.;
1610 return;
1611 }
1612
1613 frac = (wavlen-gd->wavlen[gd->cAxis][ind])/(gd->wavlen[gd->cAxis][ind+1]-gd->wavlen[gd->cAxis][ind]);
1614 ASSERT( frac > 0.-10.*DBL_EPSILON && frac < 1.+10.*DBL_EPSILON );
1615 nre = (1.-frac)*gd->n[gd->cAxis][ind].real() + frac*gd->n[gd->cAxis][ind+1].real();
1616 ASSERT( nre > 0. );
1617 nim = (1.-frac)*gd->n[gd->cAxis][ind].imag() + frac*gd->n[gd->cAxis][ind+1].imag();
1618 ASSERT( nim > 0. );
1619 nr1 = (1.-frac)*gd->nr1[gd->cAxis][ind] + frac*gd->nr1[gd->cAxis][ind+1];
1620 ASSERT( fabs(nre-1.-nr1) < 10.*max(nre,1.)*DBL_EPSILON );
1621
1622 /* size in micron, area in cm^2 */
1623 area = PI*pow2(sd->cSize)*1.e-8;
1624
1625 x = sd->cSize/wavlen*2.*PI;
1626
1627 /* note that in the following, only nre, nim are used in sinpar
1628 * and also nr1 in anomalous diffraction limit */
1629
1630 sinpar(nre,nim,x,&qext,&qphase,&qscatt,&ctbrqs,&qback,&iflag);
1631
1632 /* iflag=0 normal exit, 1 failure to converge, 2 not even tried
1633 * for exit 1,2, see whether anomalous diffraction is available */
1634
1635 if( iflag == 0 )
1636 {
1637 *error = 0;
1638 *cs_abs = area*(qext - qscatt);
1639 *cs_sct = area*qscatt;
1640 *cosb = ctbrqs/qscatt;
1641 }
1642 else
1643 {
1644 /* anomalous diffraction -- x >> 1 and |m-1| << 1 but any phase shift */
1645 if( x >= 100. && sqrt(nr1*nr1+nim*nim) <= 0.001 )
1646 {
1647 delta = -nr1;
1648 beta = nim;
1649
1650 anomal(x,&aqext,&aqabs,&aqphas,&xistar,delta,beta);
1651
1652 /* cosb is invalid */
1653 *error = 1;
1654 *cs_abs = area*aqabs;
1655 *cs_sct = area*(aqext - aqabs);
1656 *cosb = -2.;
1657 }
1658 /* nothing works */
1659 else
1660 {
1661 *error = 2;
1662 *cs_abs = -1.;
1663 *cs_sct = -1.;
1664 *cosb = -2.;
1665 }
1666 }
1667 if( *error < 2 )
1668 {
1669 if( *cs_abs <= 0. || *cs_sct <= 0. )
1670 {
1671 fprintf( ioQQQ, " illegal opacity found: wavl=%.4e micron," , wavlen );
1672 fprintf( ioQQQ, " abs_cs=%.2e, sct_cs=%.2e\n" , *cs_abs , *cs_sct );
1673 fprintf( ioQQQ, " please check refractive index file...\n" );
1675 }
1676 }
1677 if( *error < 1 )
1678 {
1679 if( fabs(*cosb) > 1.+10.*DBL_EPSILON )
1680 {
1681 fprintf( ioQQQ, " illegal asymmetry parameter found: wavl=%.4e micron," , wavlen );
1682 fprintf( ioQQQ, " cosb=%.2e\n" , *cosb );
1683 fprintf( ioQQQ, " please check refractive index file...\n" );
1685 }
1686 }
1687
1688 return;
1689}
1690
1691
1692/* this routine calculates the absorption cross sections of carbonaceous grains, it is based on Eq. 2 of:
1693 * >>refer grain physics Li, A., & Draine, B.T. 2001 ApJ, 554, 778 */
1694STATIC void ld01_fun(/*@in@*/ void(*pah_fun)(double,const sd_data*,const grain_data[],double*,double*,double*,int*),
1695 /*@in@*/ double q_gra, /* defined in LD01 */
1696 /*@in@*/ double wmin, /* below wmin use pure graphite */
1697 /*@in@*/ double wavl, /* in micron */
1698 /*@in@*/ const sd_data *sd,
1699 /*@in@*/ const grain_data gdArr[], /* gdArr[2] */
1700 /*@out@*/ double *cs_abs,
1701 /*@out@*/ double *cs_sct,
1702 /*@out@*/ double *cosb,
1703 /*@out@*/ int *error)
1704{
1705 DEBUG_ENTRY( "ld01_fun()" );
1706
1707 // this implements Eqs. 2 & 3 of LD01; it creates a gradual change from PAH-like behavior at
1708 // small radii (less than 50A) to graphite-like at large radii. The factor q_gra is non-zero
1709 // to include a small amount of "continuum" absorption even for very small grains.
1710
1711 const double a_xi = 50.e-4;
1712 double xi_PAH, cs_abs1, cs_abs2;
1713 if( wavl >= wmin )
1714 {
1715 (*pah_fun)(wavl,sd,&gdArr[0],&cs_abs1,cs_sct,cosb,error);
1716 xi_PAH = (1.-q_gra)*min(1.,pow3(a_xi/sd->cSize));
1717 }
1718 else
1719 {
1720 cs_abs1 = 0.;
1721 xi_PAH = 0.;
1722 }
1723 // ignore cs_sct, cosb, and error from pah2_fun and return the graphite ones instead.
1724 // pah2_fun never returns errors and the other two values are ignored by the upstream code
1725 mie_cs(wavl,sd,&gdArr[1],&cs_abs2,cs_sct,cosb,error);
1726 *cs_abs = xi_PAH*cs_abs1 + (1.-xi_PAH)*cs_abs2;
1727 return;
1728}
1729
1730
1731/* this routine calculates the absorption cross sections of carbonaceous grains, it creates a gradual
1732 * change from pah1_fun defined below at small grain radii to graphite-like behavior at large radii */
1733inline void car1_fun(double wavl, /* in micron */
1734 /*@in@*/ const sd_data *sd,
1735 /*@in@*/ const grain_data gdArr[], /* gdArr[2] */
1736 /*@out@*/ double *cs_abs,
1737 /*@out@*/ double *cs_sct,
1738 /*@out@*/ double *cosb,
1739 /*@out@*/ int *error)
1740{
1741 ld01_fun(pah1_fun,0.,0.,wavl,sd,gdArr,cs_abs,cs_sct,cosb,error);
1742}
1743
1744
1745/* this routine calculates the absorption cross sections of PAH molecules, it is based on:
1746 * >>refer grain physics Desert, F.-X., Boulanger, F., Puget, J. L. 1990, A&A, 237, 215
1747 *
1748 * the original version of this routine was written by Kevin Volk (University Of Calgary) */
1749
1750static const double pah1_strength[7] = { 1.4e-21,1.8e-21,1.2e-20,6.0e-21,4.0e-20,1.9e-20,1.9e-20 };
1751static const double pah1_wlBand[7] = { 3.3, 6.18, 7.7, 8.6, 11.3, 12.0, 13.3 };
1752static const double pah1_width[7] = { 0.024, 0.102, 0.24, 0.168, 0.086, 0.174, 0.174 };
1753
1754STATIC void pah1_fun(double wavl, /* in micron */
1755 /*@in@*/ const sd_data *sd,
1756 /*@in@*/ const grain_data *gd,
1757 /*@out@*/ double *cs_abs,
1758 /*@out@*/ double *cs_sct,
1759 /*@out@*/ double *cosb,
1760 /*@out@*/ int *error)
1761{
1762 long int j;
1763 double cval1,
1764 pah1_fun_v,
1765 term,
1766 term1,
1767 term2,
1768 term3,
1769 x;
1770
1771 const double p1 = 4.0e-18;
1772 const double p2 = 1.1e-18;
1773 const double p3 = 3.3e-21;
1774 const double p4 = 6.0e-21;
1775 const double p5 = 2.4e-21;
1776 const double wl1a = 5.0;
1777 const double wl1b = 7.0;
1778 const double wl1c = 9.0;
1779 const double wl1d = 9.5;
1780 const double wl2a = 11.0;
1781 const double delwl2 = 0.3;
1782 /* this is the rise interval for the second plateau */
1783 const double wl2b = wl2a + delwl2;
1784 const double wl2c = 15.0;
1785 const double wl3a = 3.2;
1786 const double wl3b = 3.57;
1787 const double wl3m = (wl3a+wl3b)/2.;
1788 const double wl3sig = 0.1476;
1789 const double x1 = 4.0;
1790 const double x2 = 5.9;
1791 const double x2a = RYD_INF/1.e4;
1792 const double x3 = 0.1;
1793
1794 /* grain volume */
1795 double vol = 4.*PI/3.*pow3(sd->cSize)*1.e-12;
1796 /* number of carbon atoms in PAH molecule */
1797 double xnc = floor(vol*gd->rho/(ATOMIC_MASS_UNIT*dense.AtomicWeight[ipCARBON]));
1798 /* number of hydrogen atoms in PAH molecule */
1799 /* >>chng 02 oct 18, use integral number of hydrogen atoms instead of fractional number */
1800 double xnh = floor(sqrt(6.*xnc));
1801 /* this is the hydrogen over carbon ratio in the PAH molecule */
1802 double xnhoc = xnh/xnc;
1803 /* ftoc3p3 is the feature to continuum ratio at 3.3 micron */
1804 double ftoc3p3 = 100.;
1805
1806 double csVal1 = 0.;
1807 double csVal2 = 0.;
1808
1809 DEBUG_ENTRY( "pah1_fun()" );
1810
1811 x = 1./wavl;
1812
1813 if( x >= x2a )
1814 {
1815 /* >>chng 02 oct 18, use atomic cross sections for energies larger than 1 Ryd */
1816 double anu_ev = x/x2a*EVRYD;
1817
1818 /* use Hartree-Slater cross sections */
1820
1821 term1 = t_ADfA::Inst().phfit(ipHYDROGEN+1,1,1,anu_ev);
1822 term2 = 0.;
1823 for( j=1; j <= 3; ++j )
1824 term2 += t_ADfA::Inst().phfit(ipCARBON+1,6,j,anu_ev);
1825
1826 csVal1 = (xnh*term1 + xnc*term2)*1.e-18;
1827 }
1828
1829 if( x <= 2.*x2a )
1830 {
1831 cval1 = log(sqrt(xnc)*ftoc3p3/1.2328)/12.2;
1832
1833 term = pow2(min(x,x1))*(3.*x1 - 2.*min(x,x1))/pow3(x1);
1834
1835 term1 = (0.1*x + 0.41)*pow2(max(x-x2,0.));
1836
1837 /* The following is an exponential cut-off in the continuum for
1838 * wavelengths larger than 2500 Angstroms; it is exponential in
1839 * wavelength so it varies as 1/x here. This replaces the
1840 * sharp cut-off at 8000 Angstroms in the original paper.
1841 *
1842 * This choice of continuum shape is also arbitrary. The continuum
1843 * is never observed at these wavelengths. For the "standard" ratio
1844 * value at 3.3 microns the continuum level in the optical is not that
1845 * much smaller than it was in the original paper. If one wants to
1846 * exactly reproduce the original optical opacity, one can change
1847 * the x1 value to a value of 1.125. Then there will be a discontinuity
1848 * in the cross-section at 8000 Angstroms.
1849 *
1850 * My judgement was that the flat cross-section in the optical used by
1851 * Desert, Boulander, and Puget (1990) is just a rough value that is not
1852 * based upon much in the way of direct observations, and so I could
1853 * change the cross-section at wavelengths above 2500 Angstroms. It is
1854 * likely that one should really build in the ERE somehow, judging from
1855 * the spectrum of the Red Rectangle, and there is no trace of this in
1856 * the original paper. The main concern in adding this exponential
1857 * drop-off in the continuum was to have a finite infrared continuum
1858 * between the features. */
1859 term2 = exp(cval1*(1. - (x1/min(x,x1))));
1860
1861 term3 = p3*exp(-pow2(x/x3))*min(x,x3)/x3;
1862
1863 csVal2 = xnc*((p1*term + p2*term1)*term2 + term3);
1864 }
1865
1866 if( x2a <= x && x <= 2.*x2a )
1867 {
1868 /* create gradual change from Desert et al to atomic cross sections */
1869 double frac = pow2(2.-x/x2a);
1870 pah1_fun_v = exp((1.-frac)*log(csVal1) + frac*log(csVal2));
1871 }
1872 else
1873 {
1874 /* one of these will be zero */
1875 pah1_fun_v = csVal1 + csVal2;
1876 }
1877
1878 /* now add in the three plateau features. the first two are based upon
1879 * >>refer grain physics Schutte, Tielens, and Allamandola (1993) ApJ, 415, 397. */
1880 if( wl1a <= wavl && wl1d >= wavl )
1881 {
1882 if( wavl < wl1b )
1883 term = p4*(wavl - wl1a)/(wl1b - wl1a);
1884 else
1885 term = ( wavl > wl1c ) ? p4*(wl1d - wavl)/(wl1d - wl1c) : p4;
1886 pah1_fun_v += term*xnc;
1887 }
1888 if( wl2a <= wavl && wl2c >= wavl )
1889 {
1890 term = ( wavl < wl2b ) ? p5*((wavl - wl2a)/delwl2) : p5*sqrt(1.-pow2((wavl-wl2a)/(wl2c-wl2a)));
1891 pah1_fun_v += term*xnc;
1892 }
1893 if( wl3m-10.*wl3sig <= wavl && wavl <= wl3m+10.*wl3sig )
1894 {
1895 /* >>chng 02 nov 08, replace top hat distribution with gaussian, PvH */
1896 term = 1.1*pah1_strength[0]*exp(-0.5*pow2((wavl-wl3m)/wl3sig));
1897 pah1_fun_v += term*xnh;
1898 }
1899
1900 /* add in the various discrete features in the infrared: 3.3, 6.2, 7.6, etc. */
1901 for( j=0; j < 7; j++ )
1902 {
1903 term1 = (wavl - pah1_wlBand[j])/pah1_width[j];
1904 term = 0.;
1905 if( j == 2 )
1906 {
1907 /* This assumes linear interpolation between the points, which are
1908 * located at -1, -0.5, +1.5, and +3 times the width, or a fine spacing that
1909 * well samples the profile. Otherwise there will be an error in the total
1910 * feature strength of order 50% */
1911 if( term1 >= -1. && term1 < -0.5 )
1912 {
1913 term = pah1_strength[j]/(3.*pah1_width[j]);
1914 term *= 2. + 2.*term1;
1915 }
1916 if( term1 >= -0.5 && term1 <= 1.5 )
1917 term = pah1_strength[j]/(3.*pah1_width[j]);
1918 if( term1 > 1.5 && term1 <= 3. )
1919 {
1920 term = pah1_strength[j]/(3.*pah1_width[j]);
1921 term *= 2. - term1*2./3.;
1922 }
1923 }
1924 else
1925 {
1926 /* This assumes linear interpolation between the points, which are
1927 * located at -2, -1, +1, and +2 times the width, or a fine spacing that
1928 * well samples the profile. Otherwise there will be an error in the total
1929 * feature strength of order 50% */
1930 if( term1 >= -2. && term1 < -1. )
1931 {
1932 term = pah1_strength[j]/(3.*pah1_width[j]);
1933 term *= 2. + term1;
1934 }
1935 if( term1 >= -1. && term1 <= 1. )
1936 term = pah1_strength[j]/(3.*pah1_width[j]);
1937 if( term1 > 1. && term1 <= 2. )
1938 {
1939 term = pah1_strength[j]/(3.*pah1_width[j]);
1940 term *= 2. - term1;
1941 }
1942 }
1943 if( j == 0 || j > 2 )
1944 term *= xnhoc;
1945 pah1_fun_v += term*xnc;
1946 }
1947
1948 *cs_abs = pah1_fun_v;
1949 /* the next two numbers are completely arbitrary, but the code requires them... */
1950 /* >>chng 02 oct 18, cs_sct was 1.e-30, but this is very high for X-ray photons, PvH */
1951 *cs_sct = 0.1*pah1_fun_v;
1952 *cosb = 0.;
1953 *error = 0;
1954
1955 return;
1956}
1957
1958
1959/* this routine calculates the absorption cross sections of carbonaceous grains, it is based on Eq. 2 of:
1960 * >>refer grain physics Li, A., & Draine, B.T. 2001 ApJ, 554, 778 */
1961inline void car2_fun(double wavl, /* in micron */
1962 /*@in@*/ const sd_data *sd,
1963 /*@in@*/ const grain_data gdArr[], /* gdArr[2] */
1964 /*@out@*/ double *cs_abs,
1965 /*@out@*/ double *cs_sct,
1966 /*@out@*/ double *cosb,
1967 /*@out@*/ int *error)
1968{
1969 ld01_fun(pah2_fun,0.01,1./17.25,wavl,sd,gdArr,cs_abs,cs_sct,cosb,error);
1970}
1971
1972
1973// these values are taken from Table 1 of LD01
1974static const double pah2_wavl[14] = { 0.0722, 0.2175, 3.3, 6.2, 7.7, 8.6, 11.3,
1975 11.9, 12.7, 16.4, 18.3, 21.2, 23.1, 26.0 };
1976static const double pah2_width[14] = { 0.195, 0.217, 0.012, 0.032, 0.091, 0.047, 0.018,
1977 0.025, 0.024, 0.010, 0.036, 0.038, 0.046, 0.69 };
1978static const double pah2n_strength[14] = { 7.97e-13, 1.23e-13, 1.97e-18, 1.96e-19, 6.09e-19, 3.47e-19, 4.27e-18,
1979 7.27e-19, 1.67e-18, 5.52e-20, 6.04e-20, 1.08e-19, 2.78e-20, 1.52e-19 };
1980static const double pah2c_strength[14] = { 7.97e-13, 1.23e-13, 4.47e-19, 1.57e-18, 5.48e-18, 2.42e-18, 4.00e-18,
1981 6.14e-19, 1.49e-18, 5.52e-20, 6.04e-20, 1.08e-19, 2.78e-20, 1.52e-19 };
1982
1983/* this routine calculates the absorption cross sections of PAH molecules, it is based on Eqs. 6-11 of:
1984 * >>refer grain physics Li, A., & Draine, B.T. 2001 ApJ, 554, 778 */
1985STATIC void pah2_fun(double wavl, /* in micron */
1986 /*@in@*/ const sd_data *sd,
1987 /*@in@*/ const grain_data *gd,
1988 /*@out@*/ double *cs_abs,
1989 /*@out@*/ double *cs_sct,
1990 /*@out@*/ double *cosb,
1991 /*@out@*/ int *error)
1992{
1993 DEBUG_ENTRY( "pah2_fun()" );
1994
1995 // these choices give the best fit to the observed spectrum of the diffuse ISM
1996 // setting all these to 1 reproduces the laboratory measured band strengths
1997 const double E62 = 3.;
1998 const double E77 = 2.;
1999 const double E86 = 2.;
2000
2001 // grain volume
2002 double vol = 4.*PI/3.*pow3(sd->cSize)*1.e-12;
2003 // number of carbon atoms in PAH molecule
2004 double xnc = vol*gd->rho/(ATOMIC_MASS_UNIT*dense.AtomicWeight[ipCARBON]);
2005 // this is the hydrogen over carbon ratio in the PAH molecule
2006 double xnhoc;
2007 // this is Eq. 4 of LD01
2008 if( xnc <= 25. )
2009 xnhoc = 0.5;
2010 else if( xnc <= 100. )
2011 xnhoc = 2.5/sqrt(xnc);
2012 else
2013 xnhoc = 0.25;
2014
2015 double x = 1./wavl;
2016
2017 double pah2_fun_v = 0.;
2018 if( x < 3.3 )
2019 {
2020 double M = ( xnc <= 40. ) ? 0.3*xnc : 0.4*xnc;
2021 // calculate cutoff wavelength in micron, this is Eq. A3 of LD01
2022 double cutoff;
2023 if( gd->charge == 0 )
2024 cutoff = 1./(3.804/sqrt(M) + 1.052);
2025 else
2026 cutoff = 1./(2.282/sqrt(M) + 0.889);
2027 double y = cutoff/wavl;
2028 // this is Eq. A2 of LD01
2029 double cutoff_fun = atan(1.e3*pow3(y-1.)/y)/PI + 0.5;
2030 // this is Eq. 11 of LD01
2031 pah2_fun_v = 34.58*pow( 10., -18.-3.431/x )*cutoff_fun;
2032 for( int j=2; j < 14; ++j )
2033 {
2034 double strength = ( gd->charge == 0 ) ? pah2n_strength[j] : pah2c_strength[j];
2035 if( j == 2 )
2036 strength *= xnhoc;
2037 else if( j == 3 )
2038 strength *= E62;
2039 else if( j == 4 )
2040 strength *= E77;
2041 else if( j == 5 )
2042 strength *= E86*xnhoc;
2043 else if( j == 6 || j == 7 || j == 8 )
2044 strength *= xnhoc/3.;
2045 pah2_fun_v += Drude( wavl, pah2_wavl[j], pah2_width[j], strength );
2046 }
2047 }
2048 else if( x < 5.9 )
2049 {
2050 // this is Eq. 10 of LD01, strength is identical for charged and neutral PAHs
2051 pah2_fun_v = Drude( wavl, pah2_wavl[1], pah2_width[1], pah2n_strength[1] );
2052 pah2_fun_v += (1.8687 + 0.1905*x)*1.e-18;
2053 }
2054 else if( x < 7.7 )
2055 {
2056 // this is Eq. 9 of LD01, strength is identical for charged and neutral PAHs
2057 pah2_fun_v = Drude( wavl, pah2_wavl[1], pah2_width[1], pah2n_strength[1] );
2058 double y = x - 5.9;
2059 pah2_fun_v += (1.8687 + 0.1905*x + pow2(y)*(0.4175 + 0.04370*y))*1.e-18;
2060 }
2061 else if( x < 10. )
2062 {
2063 // this is Eq. 8 of LD01
2064 pah2_fun_v = (((-0.1057*x + 2.950)*x - 24.367)*x + 66.302)*1.e-18;
2065 }
2066 else if( x < 15. )
2067 {
2068 // this is Eq. 7 of LD01, strength is identical for charged and neutral PAHs
2069 pah2_fun_v = Drude( wavl, pah2_wavl[0], pah2_width[0], pah2n_strength[0] );
2070 pah2_fun_v += (-3. + 1.35*x)*1.e-18;
2071 }
2072 else if( x < 17.26 )
2073 {
2074 // this is Eq. 6 of LD01
2075 pah2_fun_v = (126.0 - 6.4943*x)*1.e-18;
2076 }
2077 else
2078 // this routine should never be called for wavelengths shorter than 1/17.25 micron
2079 // graphite opacities should be used in that case; this is handled in car2_fun()
2080 TotalInsanity();
2081
2082 // normalize cross section per PAH molecule
2083 pah2_fun_v *= xnc;
2084
2085 *cs_abs = pah2_fun_v;
2086 // the next two numbers are completely arbitrary
2087 *cs_sct = 0.1*pah2_fun_v;
2088 *cosb = 0.;
2089 *error = 0;
2090
2091 return;
2092}
2093
2094
2095/* this routine calculates the absorption cross sections of carbonaceous grains, it is based on Eqs. 5-7 of:
2096 * >>refer grain physics Draine, B.T., & Li, A., 2007 ApJ, 657, 810 */
2097inline void car3_fun(double wavl, /* in micron */
2098 /*@in@*/ const sd_data *sd,
2099 /*@in@*/ const grain_data gdArr[], /* gdArr[2] */
2100 /*@out@*/ double *cs_abs,
2101 /*@out@*/ double *cs_sct,
2102 /*@out@*/ double *cosb,
2103 /*@out@*/ int *error)
2104{
2105 ld01_fun(pah3_fun,0.01,1./17.25,wavl,sd,gdArr,cs_abs,cs_sct,cosb,error);
2106}
2107
2108
2109// these values are taken from Table 1 of DL07
2110static const double pah3_wavl[30] = { 0.0722, 0.2175, 1.05, 1.26, 1.905, 3.3, 5.27, 5.7, 6.22, 6.69, 7.417, 7.598,
2111 7.85, 8.33, 8.61, 10.68, 11.23, 11.33, 11.99, 12.62, 12.69, 13.48, 14.19,
2112 15.9, 16.45, 17.04, 17.375, 17.87, 18.92, 15. };
2113static const double pah3_width[30] = { 0.195, 0.217, 0.055, 0.11, 0.09, 0.012, 0.034, 0.035, 0.03, 0.07, 0.126,
2114 0.044, 0.053, 0.052, 0.039, 0.02, 0.012, 0.032, 0.045, 0.042, 0.013, 0.04,
2115 0.025, 0.02, 0.014, 0.065, 0.012, 0.016, 0.1, 0.8 };
2116static const double pah3n_strength[30] = { 7.97e-13, 1.23e-13, 0., 0., 0., 3.94e-18, 2.5e-20, 4.e-20, 2.94e-19,
2117 7.35e-20, 2.08e-19, 1.81e-19, 2.19e-19, 6.94e-20, 2.78e-19, 3.e-21,
2118 1.89e-19, 5.2e-19, 2.42e-19, 3.5e-19, 1.3e-20, 8.e-20, 4.5e-21,
2119 4.e-22, 5.e-21, 2.22e-20, 1.1e-21, 6.7e-22, 1.e-21, 5.e-19 };
2120static const double pah3c_strength[30] = { 7.97e-13, 1.23e-13, 2.e-16, 7.8e-17, -1.465e-18, 8.94e-19, 2.e-19,
2121 3.2e-19, 2.35e-18, 5.9e-19, 1.81e-18, 1.63e-18, 1.97e-18, 4.8e-19,
2122 1.94e-18, 3.e-21, 1.77e-19, 4.9e-19, 2.05e-19, 3.1e-19, 1.3e-20,
2123 8.e-20, 4.5e-21, 4.e-22, 5.e-21, 2.22e-20, 1.1e-21, 6.7e-22, 1.7e-21,
2124 5.e-19 };
2125// should we multiply the strength with H/C ?
2126static const bool pah3_hoc[30] = { false, false, false, false, false, true, false, false, false, false, false,
2127 false, false, true, true, true, true, true, true, true, true, true, false,
2128 false, false, false, false, false, false, false };
2129
2130/* this routine calculates the absorption cross sections of PAH molecules, it is based on
2131 * >>refer grain physics Draine, B.T., & Li, A., 2007 ApJ, 657, 810 */
2132STATIC void pah3_fun(double wavl, /* in micron */
2133 /*@in@*/ const sd_data *sd,
2134 /*@in@*/ const grain_data *gd,
2135 /*@out@*/ double *cs_abs,
2136 /*@out@*/ double *cs_sct,
2137 /*@out@*/ double *cosb,
2138 /*@out@*/ int *error)
2139{
2140 DEBUG_ENTRY( "pah3_fun()" );
2141
2142 // grain volume
2143 double vol = 4.*PI/3.*pow3(sd->cSize)*1.e-12;
2144 // number of carbon atoms in PAH molecule
2145 double xnc = vol*gd->rho/(ATOMIC_MASS_UNIT*dense.AtomicWeight[ipCARBON]);
2146 // this is the hydrogen over carbon ratio in the PAH molecule
2147 double xnhoc;
2148 // this is Eq. 4 of DL07
2149 if( xnc <= 25. )
2150 xnhoc = 0.5;
2151 else if( xnc <= 100. )
2152 xnhoc = 2.5/sqrt(xnc);
2153 else
2154 xnhoc = 0.25;
2155
2156 double x = 1./wavl;
2157
2158 // this is Eq. 2 of DL07
2159 double pah3_fun_v = ( gd->charge == 0 ) ? 0. : 3.5*pow(10.,-19.-1.45/x)*exp(-0.1*pow2(x));
2160
2161 if( x < 3.3 )
2162 {
2163 double M = ( xnc <= 40. ) ? 0.3*xnc : 0.4*xnc;
2164 // calculate cutoff wavelength in micron, this is Eq. A3 of LD01
2165 double cutoff;
2166 if( gd->charge == 0 )
2167 cutoff = 1./(3.804/sqrt(M) + 1.052);
2168 else
2169 cutoff = 1./(2.282/sqrt(M) + 0.889);
2170 double y = cutoff/wavl;
2171 // this is Eq. A2 of LD01
2172 double cutoff_fun = atan(1.e3*pow3(y-1.)/y)/PI + 0.5;
2173 // this is Eq. 11 of LD01
2174 pah3_fun_v += 34.58*pow( 10., -18.-3.431/x )*cutoff_fun;
2175 for( int j=2; j < 30; ++j )
2176 {
2177 double strength = ( gd->charge == 0 ) ? pah3n_strength[j] : pah3c_strength[j];
2178 if( pah3_hoc[j] )
2179 strength *= xnhoc;
2180 pah3_fun_v += Drude( wavl, pah3_wavl[j], pah3_width[j], strength );
2181 }
2182 }
2183 else if( x < 5.9 )
2184 {
2185 // this is Eq. 10 of LD01, strength is identical for charged and neutral PAHs
2186 pah3_fun_v += Drude( wavl, pah3_wavl[1], pah3_width[1], pah3n_strength[1] );
2187 pah3_fun_v += (1.8687 + 0.1905*x)*1.e-18;
2188 }
2189 else if( x < 7.7 )
2190 {
2191 // this is Eq. 9 of LD01, strength is identical for charged and neutral PAHs
2192 pah3_fun_v += Drude( wavl, pah3_wavl[1], pah3_width[1], pah3n_strength[1] );
2193 double y = x - 5.9;
2194 pah3_fun_v += (1.8687 + 0.1905*x + pow2(y)*(0.4175 + 0.04370*y))*1.e-18;
2195 }
2196 else if( x < 10. )
2197 {
2198 // this is Eq. 8 of LD01
2199 pah3_fun_v += (((-0.1057*x + 2.950)*x - 24.367)*x + 66.302)*1.e-18;
2200 }
2201 else if( x < 15. )
2202 {
2203 // this is Eq. 7 of LD01, strength is identical for charged and neutral PAHs
2204 pah3_fun_v += Drude( wavl, pah3_wavl[0], pah3_width[0], pah3n_strength[0] );
2205 pah3_fun_v += (-3. + 1.35*x)*1.e-18;
2206 }
2207 else if( x < 17.26 )
2208 {
2209 // this is Eq. 6 of LD01
2210 pah3_fun_v += (126.0 - 6.4943*x)*1.e-18;
2211 }
2212 else
2213 // this routine should never be called for wavelengths shorter than 1/17.25 micron
2214 // graphite opacities should be used in that case; this is handled in car3_fun()
2215 TotalInsanity();
2216
2217 // normalize cross section per PAH molecule
2218 pah3_fun_v *= xnc;
2219
2220 *cs_abs = pah3_fun_v;
2221 // the next two numbers are completely arbitrary
2222 *cs_sct = 0.1*pah3_fun_v;
2223 *cosb = 0.;
2224 *error = 0;
2225
2226 return;
2227}
2228
2229
2230// Drude: helper function to calculate Drude profile, return value is cross section in cm^2/C
2231inline double Drude(double lambda, // wavelength (in micron)
2232 double lambdac, // central wavelength of feature (in micron)
2233 double gamma, // width of the feature (dimensionless)
2234 double sigma) // strength of the feature (in cm/C)
2235{
2236 double x = lambda/lambdac;
2237 // this is Eq. 12 of LD01
2238 return 2.e-4/PI*gamma*lambdac*sigma/(pow2(x - 1./x) + pow2(gamma));
2239}
2240
2241
2242STATIC void tbl_fun(double wavl, /* in micron */
2243 /*@in@*/ const sd_data *sd, /* NOT USED -- MUST BE KEPT FOR COMPATIBILITY */
2244 /*@in@*/ const grain_data *gd,
2245 /*@out@*/ double *cs_abs,
2246 /*@out@*/ double *cs_sct,
2247 /*@out@*/ double *cosb,
2248 /*@out@*/ int *error)
2249{
2250 bool lgOutOfBounds;
2251 long int ind;
2252 double anu = WAVNRYD/wavl*1.e4;
2253
2254 DEBUG_ENTRY( "tbl_fun()" );
2255
2256 /* >>chng 02 nov 17, add this test to prevent warning that this var not used */
2257 if( sd == NULL )
2258 TotalInsanity();
2259
2261
2262 find_arr(anu,gd->opcAnu,gd->nOpcData,&ind,&lgOutOfBounds);
2263 if( !lgOutOfBounds )
2264 {
2265 double a1g;
2266 double frac = log(anu/gd->opcAnu[ind])/log(gd->opcAnu[ind+1]/gd->opcAnu[ind]);
2267
2268 *cs_abs = exp((1.-frac)*log(gd->opcData[0][ind])+frac*log(gd->opcData[0][ind+1]));
2269 ASSERT( *cs_abs > 0. );
2270 if( gd->nOpcCols > 1 )
2271 *cs_sct = exp((1.-frac)*log(gd->opcData[1][ind])+frac*log(gd->opcData[1][ind+1]));
2272 else
2273 *cs_sct = 0.1*(*cs_abs);
2274 ASSERT( *cs_sct > 0. );
2275 if( gd->nOpcCols > 2 )
2276 a1g = exp((1.-frac)*log(gd->opcData[2][ind])+frac*log(gd->opcData[2][ind+1]));
2277 else
2278 a1g = 1.;
2279 ASSERT( a1g > 0. );
2280 *cosb = 1. - a1g;
2281 *error = 0;
2282 }
2283 else
2284 {
2285 *cs_abs = -1.;
2286 *cs_sct = -1.;
2287 *cosb = -2.;
2288 *error = 3;
2289 }
2290 return;
2291}
2292
2293
2294STATIC double size_distr(double size,
2295 /*@in@*/ const sd_data *sd)
2296{
2297 bool lgOutOfBounds;
2298 long ind;
2299 double frac,
2300 res,
2301 x;
2302
2303 DEBUG_ENTRY( "size_distr()" );
2304
2305 if( size >= sd->lim[ipBLo] && size <= sd->lim[ipBHi] )
2306 switch( sd->sdCase )
2307 {
2308 case SD_SINGLE_SIZE:
2309 case SD_NR_CARBON:
2310 res = 1.; /* should really not be used in this case */
2311 break;
2312 case SD_POWERLAW:
2313 /* simple powerlaw */
2314 case SD_EXP_CUTOFF1:
2315 case SD_EXP_CUTOFF2:
2316 case SD_EXP_CUTOFF3:
2317 /* powerlaw with exponential cutoff, inspired by Greenberg (1978)
2318 * Cosmic Dust, ed. J.A.M. McDonnell, Wiley, p. 187 */
2319 res = pow(size,sd->a[ipExp]);
2320 if( sd->a[ipBeta] < 0. )
2321 res /= (1. - sd->a[ipBeta]*size);
2322 else if( sd->a[ipBeta] > 0. )
2323 res *= (1. + sd->a[ipBeta]*size);
2324 if( size < sd->a[ipBLo] && sd->a[ipSLo] > 0. )
2325 res *= exp(-powi((sd->a[ipBLo]-size)/sd->a[ipSLo],nint(sd->a[ipAlpha])));
2326 if( size > sd->a[ipBHi] && sd->a[ipSHi] > 0. )
2327 res *= exp(-powi((size-sd->a[ipBHi])/sd->a[ipSHi],nint(sd->a[ipAlpha])));
2328 break;
2329 case SD_LOG_NORMAL:
2330 x = log(size/sd->a[ipGCen])/sd->a[ipGSig];
2331 res = exp(-0.5*pow2(x))/size;
2332 break;
2333 case SD_LIN_NORMAL:
2334 x = (size-sd->a[ipGCen])/sd->a[ipGSig];
2335 res = exp(-0.5*pow2(x))/size;
2336 break;
2337 case SD_TABLE:
2338 find_arr(log(size),sd->ln_a,sd->npts,&ind,&lgOutOfBounds);
2339 if( lgOutOfBounds )
2340 {
2341 fprintf( ioQQQ, " size distribution table has insufficient range\n" );
2342 fprintf( ioQQQ, " requested size: %.5f table range %.5f - %.5f\n",
2343 size, exp(sd->ln_a[0]), exp(sd->ln_a[sd->npts-1]) );
2345 }
2346 frac = (log(size)-sd->ln_a[ind])/(sd->ln_a[ind+1]-sd->ln_a[ind]);
2347 ASSERT( frac > 0.-10.*DBL_EPSILON && frac < 1.+10.*DBL_EPSILON );
2348 res = (1.-frac)*sd->ln_a4dNda[ind] + frac*sd->ln_a4dNda[ind+1];
2349 /* convert from a^4*dN/da to dN/da */
2350 res = exp(res)/POW4(size);
2351 break;
2352 case SD_ILLEGAL:
2353 default:
2354 fprintf( ioQQQ, " insane case for grain size distribution: %d\n" , sd->sdCase );
2355 ShowMe();
2357 }
2358 else
2359 res = 0.;
2360 return res;
2361}
2362
2363
2364/* search for upper/lower limit lim of size distribution such that
2365 * lim^4 * dN/da(lim) < rel_cutoff * ref^4 * dN/da(ref)
2366 * the initial estimate of lim = ref + step
2367 * step may be both positive (upper limit) or negative (lower limit) */
2368STATIC double search_limit(double ref,
2369 double step,
2370 double rel_cutoff,
2371 // do NOT use sd_data* since that would trash sd.lim[ipBLo]
2372 sd_data sd)
2373{
2374 long i;
2375 double f1,
2376 f2,
2377 fmid,
2378 renorm,
2379 x1,
2380 x2 = DBL_MAX,
2381 xmid = DBL_MAX;
2382
2383 /* TOLER is the relative accuracy with which lim is determined */
2384 const double TOLER = 1.e-6;
2385
2386 DEBUG_ENTRY( "search_limit()" );
2387
2388 /* sanity check */
2389 ASSERT( rel_cutoff > 0. && rel_cutoff < 1. );
2390
2391 if( step == 0. )
2392 {
2393 return ref;
2394 }
2395
2396 /* these need to be set in order for size_distr to work...
2397 * NB - since this is a local copy of sd, it will not
2398 * upset anything in the calling routine */
2399 sd.lim[ipBLo] = 0.;
2400 sd.lim[ipBHi] = DBL_MAX;
2401
2402 x1 = ref;
2403 /* previous assert guarantees that f1 > 0. */
2404 f1 = -log(rel_cutoff);
2405 renorm = f1 - log(POW4(x1)*size_distr(x1,&sd));
2406
2407 /* bracket solution */
2408 f2 = 1.;
2409 for( i=0; i < 20 && f2 > 0.; ++i )
2410 {
2411 x2 = max(ref + step,SMALLEST_GRAIN);
2412 f2 = log(POW4(x2)*size_distr(x2,&sd)) + renorm;
2413 if( f2 >= 0. )
2414 {
2415 x1 = x2;
2416 f1 = f2;
2417 }
2418 step *= 2.;
2419 }
2420 if( f2 > 0. )
2421 {
2422 fprintf( ioQQQ, " Could not bracket solution\n" );
2424 }
2425
2426 /* do bisection search */
2427 while( 2.*fabs(x1-x2)/(x1+x2) > TOLER )
2428 {
2429 xmid = (x1+x2)/2.;
2430 fmid = log(POW4(xmid)*size_distr(xmid,&sd)) + renorm;
2431
2432 if( fmid == 0. )
2433 break;
2434
2435 if( f1*fmid > 0. )
2436 {
2437 x1 = xmid;
2438 f1 = fmid;
2439 }
2440 else
2441 {
2442 x2 = xmid;
2443 f2 = fmid;
2444 }
2445 }
2446 return (x1+x2)/2.;
2447}
2448
2449/* calculate the inverse attenuation length for given refractive index data */
2450STATIC void mie_calc_ial(/*@in@*/ const grain_data *gd,
2451 long int n,
2452 /*@out@*/ vector<double>& invlen, /* invlen[n] */
2453 /*@in@*/ const char *chString,
2454 /*@in@*/ bool *lgWarning)
2455{
2456 bool lgErrorOccurred=true,
2457 lgOutOfBounds;
2458 long int i,
2459 ind,
2460 j;
2461 double frac,
2462 InvDep,
2463 nim,
2464 wavlen;
2465
2466 DEBUG_ENTRY( "mie_calc_ial()" );
2467
2468 /* sanity check */
2469 ASSERT( gd->rfiType == RFI_TABLE );
2470
2471 vector<int> ErrorIndex( rfield.nupper );
2472
2473 for( i=0; i < n; i++ )
2474 {
2475 wavlen = WAVNRYD/rfield.anu[i]*1.e4;
2476
2477 ErrorIndex[i] = 0;
2478 lgErrorOccurred = false;
2479 invlen[i] = 0.;
2480
2481 for( j=0; j < gd->nAxes; j++ )
2482 {
2483 /* first interpolate optical constants */
2484 find_arr(wavlen,gd->wavlen[j],gd->ndata[j],&ind,&lgOutOfBounds);
2485 if( lgOutOfBounds )
2486 {
2487 ErrorIndex[i] = 3;
2488 lgErrorOccurred = true;
2489 invlen[i] = 0.;
2490 break;
2491 }
2492 frac = (wavlen-gd->wavlen[j][ind])/(gd->wavlen[j][ind+1]-gd->wavlen[j][ind]);
2493 nim = (1.-frac)*gd->n[j][ind].imag() + frac*gd->n[j][ind+1].imag();
2494 /* this is the inverse of the photon attenuation depth,
2495 * >>refer grain physics Weingartner & Draine, 2000, ApJ, ... */
2496 InvDep = PI4*nim/wavlen*1.e4;
2497 ASSERT( InvDep > 0. );
2498
2499 invlen[i] += InvDep*gd->wt[j];
2500 }
2501 }
2502
2503 if( lgErrorOccurred )
2504 {
2505 mie_repair(chString,n,3,3,rfield.anu,&invlen[0],ErrorIndex,false,lgWarning);
2506 }
2507
2508 return;
2509}
2510
2511
2512/* this is the number of x-values we use for extrapolating functions */
2513const int NPTS_DERIV = 8;
2515
2516/* extrapolate/interpolate mie data to fill in the blanks */
2517STATIC void mie_repair(/*@in@*/ const char *chString,
2518 long int n,
2519 int val,
2520 int del,
2521 /*@in@*/ const double anu[], /* anu[n] */
2522 double data[], /* data[n] */
2523 /*@in@*/ vector<int>& ErrorIndex, /* ErrorIndex[n] */
2524 bool lgRound,
2525 /*@in@*/ bool *lgWarning)
2526{
2527 bool lgExtrapolate,
2528 lgVerbose;
2529 long int i1,
2530 i2,
2531 ind1,
2532 ind2,
2533 j;
2534 double dx,
2535 sgn,
2536 slp1,
2537 xlg1,
2538 xlg2,
2539 y1lg1,
2540 y1lg2;
2541
2542 /* interpolating over more that this number of points results in a warning */
2543 const long BIG_INTERPOLATION = 10;
2544
2545 DEBUG_ENTRY( "mie_repair()" );
2546
2547 lgVerbose = ( chString[0] != '\0' );
2548
2549 for( ind1=0; ind1 < n; )
2550 {
2551 if( ErrorIndex[ind1] == val )
2552 {
2553 /* search for region with identical error index */
2554 ind2 = ind1;
2555 while( ind2 < n && ErrorIndex[ind2] == val )
2556 ind2++;
2557
2558 if( lgVerbose )
2559 fprintf( ioQQQ, " %s", chString );
2560
2561 if( ind1 == 0 )
2562 {
2563 /* low energy extrapolation */
2564 i1 = ind2;
2565 i2 = ind2+NPTS_DERIV-1;
2566 lgExtrapolate = true;
2567 sgn = +1.;
2568 if( lgVerbose )
2569 {
2570 fprintf( ioQQQ, " extrapolated below %.4e Ryd\n",anu[i1] );
2571 }
2572 }
2573 else if( ind2 == n )
2574 {
2575 /* high energy extrapolation */
2576 i1 = ind1-NPTS_DERIV;
2577 i2 = ind1-1;
2578 lgExtrapolate = true;
2579 sgn = -1.;
2580 if( lgVerbose )
2581 {
2582 fprintf( ioQQQ, " extrapolated above %.4e Ryd\n",anu[i2] );
2583 }
2584 }
2585 else
2586 {
2587 /* interpolation */
2588 i1 = ind1-1;
2589 i2 = ind2;
2590 lgExtrapolate = false;
2591 sgn = 0.;
2592 if( lgVerbose )
2593 {
2594 fprintf( ioQQQ, " interpolated between %.4e and %.4e Ryd\n",
2595 anu[i1],anu[i2] );
2596 }
2597 if( i2-i1-1 > BIG_INTERPOLATION )
2598 {
2599 if( lgVerbose )
2600 {
2601 fprintf( ioQQQ, " ***Warning: extensive interpolation used\n" );
2602 }
2603 *lgWarning = true;
2604 }
2605 }
2606
2607 if( i1 < 0 || i2 >= n )
2608 {
2609 fprintf( ioQQQ, " Insufficient data for extrapolation\n" );
2611 }
2612
2613 xlg1 = log(anu[i1]);
2614 y1lg1 = log(data[i1]);
2615 /* >>chng 01 jul 30, replace simple-minded extrapolation with more robust routine, PvH */
2616 if( lgExtrapolate )
2617 slp1 = mie_find_slope(anu,data,ErrorIndex,i1,i2,val,lgVerbose,lgWarning);
2618 else
2619 {
2620 xlg2 = log(anu[i2]);
2621 y1lg2 = log(data[i2]);
2622 slp1 = (y1lg2-y1lg1)/(xlg2-xlg1);
2623 }
2624 if( lgRound && lgExtrapolate && sgn > 0. )
2625 {
2626 /* in low energy extrapolation, 1-g is very close to 1 and almost constant
2627 * hence slp1 is very close to 0 and can even be slightly negative
2628 * to prevent 1-g becoming greater than 1, the following is necessary */
2629 slp1 = max(slp1,0.);
2630 }
2631 /* >>chng 02 oct 22, changed from sgn*slp1 <= 0. to accomodate grey grains */
2632 else if( lgExtrapolate && sgn*slp1 < 0. )
2633 {
2634 fprintf( ioQQQ, " Unphysical value for slope in extrapolation %.6e\n", slp1 );
2635 fprintf( ioQQQ, " The most likely cause is that your refractive index or "
2636 "opacity data do not extend to low or high enough frequencies. "
2637 "See Hazy 1 for more details.\n" );
2639 }
2640
2641 for( j=ind1; j < ind2; j++ )
2642 {
2643 dx = log(anu[j]) - xlg1;
2644 data[j] = exp(y1lg1 + dx*slp1);
2645 ErrorIndex[j] -= del;
2646 }
2647
2648 ind1 = ind2;
2649 }
2650 else
2651 {
2652 ind1++;
2653 }
2654 }
2655 /* sanity check */
2656 for( j=0; j < n; j++ )
2657 {
2658 if( ErrorIndex[j] > val-del )
2659 {
2660 fprintf( ioQQQ, " Internal error in mie_repair, index=%ld, val=%d\n",j,ErrorIndex[j] );
2661 ShowMe();
2663 }
2664 }
2665 return;
2666}
2667
2668
2669STATIC double mie_find_slope(/*@in@*/ const double anu[],
2670 /*@in@*/ const double data[],
2671 /*@in@*/ const vector<int>& ErrorIndex,
2672 long i1,
2673 long i2,
2674 int val,
2675 bool lgVerbose,
2676 /*@in@*/ bool *lgWarning)
2677{
2678 long i,
2679 j,
2680 k;
2681 double s1,
2682 s2,
2683 slope,
2684 slp1[NPTS_COMB],
2685 stdev;
2686
2687 /* threshold for standard deviation in the logarithmic derivative to generate warning,
2688 * this corresponds to an uncertainty of a factor 10 for a typical extrapolation */
2689 const double LARGE_STDEV = 0.2;
2690
2691 DEBUG_ENTRY( "mie_find_slope()" );
2692
2693 /* sanity check */
2694 ASSERT( i2-i1 == NPTS_DERIV-1 );
2695 for( i=i1; i <= i2; i++ )
2696 {
2697 ASSERT( ErrorIndex[i] < val );
2698 ASSERT( anu[i] > 0. && data[i] > 0. );
2699 }
2700
2701 /* this initialization is to keep lint happy, the next statement will do the real initialization */
2702 for( i=0; i < NPTS_COMB; i++ )
2703 {
2704 slp1[i] = -DBL_MAX;
2705 }
2706
2707 k = 0;
2708 /* calculate the logarithmic derivative for every possible combination of data points */
2709 /* this loop is guaranteed to initialize all members of slp1, the first ASSERT checks this */
2710 for( i=i1; i < i2; i++ )
2711 {
2712 for( j=i+1; j <= i2; j++ )
2713 {
2714 slp1[k++] = log(data[j]/data[i])/log(anu[j]/anu[i]);
2715 }
2716 }
2717 /* sort the values; we want the median -> values for i > NPTS_COMB/2 are irrelevant */
2718 for( i=0; i <= NPTS_COMB/2; i++ )
2719 {
2720 for( j=i+1; j < NPTS_COMB; j++ )
2721 {
2722 if( slp1[i] > slp1[j] )
2723 {
2724 double xxx = slp1[i];
2725 slp1[i] = slp1[j];
2726 slp1[j] = xxx;
2727 }
2728 }
2729 }
2730 /* now calculate the median value */
2731 slope = ( NPTS_COMB%2 == 1 ) ? slp1[NPTS_COMB/2] : (slp1[NPTS_COMB/2-1] + slp1[NPTS_COMB/2])/2.;
2732
2733 /* and finally calculate the standard deviation of all slopes */
2734 s1 = s2 = 0.;
2735 for( i=0; i < NPTS_COMB; i++ )
2736 {
2737 s1 += slp1[i];
2738 s2 += pow2(slp1[i]);
2739 }
2740 /* >>chng 06 jul 12, protect against roundoff error, PvH */
2741 stdev = max(s2/(double)NPTS_COMB - pow2(s1/(double)NPTS_COMB),0.);
2742 stdev = sqrt(stdev);
2743
2744 /* print warning if standard deviation is large */
2745 if( stdev > LARGE_STDEV )
2746 {
2747 if( lgVerbose )
2748 {
2749 fprintf( ioQQQ, " ***Warning: slope for extrapolation may be unreliable\n" );
2750 }
2751 *lgWarning = true;
2752 }
2753 return slope;
2754}
2755
2756
2757/* read in the file with optical constants and other relevant information */
2758STATIC void mie_read_rfi(/*@in@*/ const char *chFile,
2759 /*@out@*/ grain_data *gd)
2760{
2761 bool lgLogData=false;
2762 long int dl,
2763 help,
2764 i,
2765 nelem,
2766 j,
2767 nridf,
2768 sgn = 0;
2769 double eps1,
2770 eps2,
2771 LargestLog,
2772 molw,
2773 nAtoms,
2774 nr,
2775 ni,
2776 tmp1,
2777 tmp2,
2778 total = 0.;
2779 char chLine[FILENAME_PATH_LENGTH_2],
2780 chWord[FILENAME_PATH_LENGTH_2];
2781 FILE *io2;
2782
2783 DEBUG_ENTRY( "mie_read_rfi()" );
2784
2785 io2 = open_data( chFile, "r", AS_LOCAL_ONLY );
2786
2787 dl = 0; /* line counter for input file */
2788
2789 /* first read magic number */
2790 mie_next_data(chFile,io2,chLine,&dl);
2791 mie_read_long(chFile,chLine,&gd->magic,true,dl);
2792 if( gd->magic != MAGIC_RFI )
2793 {
2794 fprintf( ioQQQ, " Refractive index file %s has obsolete magic number\n",chFile );
2795 fprintf( ioQQQ, " I found %ld, but expected %ld on line #%ld\n",gd->magic,MAGIC_RFI,dl );
2796 fprintf( ioQQQ, " Please replace this file with an up to date version\n" );
2798 }
2799
2800 /* get chemical formula of the grain, e.g., Mg0.4Fe0.6SiO3 */
2801 mie_next_data(chFile,io2,chLine,&dl);
2802 mie_read_word(chLine,chWord,FILENAME_PATH_LENGTH_2,false);
2803 mie_read_form(chWord,gd->elmAbun,&nAtoms,&molw);
2804
2805 /* molecular weight, in atomic units */
2806 gd->mol_weight = molw;
2807 gd->atom_weight = gd->mol_weight/nAtoms;
2808
2809 /* determine abundance of grain molecule assuming max depletion */
2810 mie_next_data(chFile,io2,chLine,&dl);
2811 mie_read_double(chFile,chLine,&gd->abun,true,dl);
2812
2813 /* default depletion of grain molecule */
2814 mie_next_data(chFile,io2,chLine,&dl);
2815 mie_read_double(chFile,chLine,&gd->depl,true,dl);
2816 if( gd->depl > 1. )
2817 {
2818 fprintf( ioQQQ, " Illegal value for default depletion in %s\n",chFile );
2819 fprintf( ioQQQ, " Line #%ld, depl=%14.6e\n",dl,gd->depl);
2821 }
2822
2823 for( nelem=0; nelem < LIMELM; nelem++ )
2824 gd->elmAbun[nelem] *= gd->abun*gd->depl;
2825
2826 /* material density, to get cross section per unit mass: rho in cgs */
2827 mie_next_data(chFile,io2,chLine,&dl);
2828 mie_read_double(chFile,chLine,&gd->rho,false,dl);
2829
2830 /* material type, determines enthalpy function: 1 -- carbonaceous, 2 -- silicate */
2831 mie_next_data(chFile,io2,chLine,&dl);
2832 mie_read_long(chFile,chLine,&help,true,dl);
2833 gd->matType = (mat_type)help;
2834 if( gd->matType >= MAT_TOP )
2835 {
2836 fprintf( ioQQQ, " Illegal value for material type in %s\n",chFile );
2837 fprintf( ioQQQ, " Line #%ld, type=%d\n",dl,gd->matType);
2839 }
2840
2841 /* work function, in Rydberg */
2842 mie_next_data(chFile,io2,chLine,&dl);
2843 mie_read_double(chFile,chLine,&gd->work,true,dl);
2844
2845 /* bandgap, in Rydberg */
2846 mie_next_data(chFile,io2,chLine,&dl);
2847 mie_read_double(chFile,chLine,&gd->bandgap,false,dl);
2848 if( gd->bandgap >= gd->work )
2849 {
2850 fprintf( ioQQQ, " Illegal value for bandgap in %s\n",chFile );
2851 fprintf( ioQQQ, " Line #%ld, bandgap=%.4e, work function=%.4e\n",dl,gd->bandgap,gd->work);
2852 fprintf( ioQQQ, " Bandgap should always be less than work function\n");
2854 }
2855
2856 /* efficiency of thermionic emission, between 0 and 1 */
2857 mie_next_data(chFile,io2,chLine,&dl);
2858 mie_read_double(chFile,chLine,&gd->therm_eff,true,dl);
2859 if( gd->therm_eff > 1.f )
2860 {
2861 fprintf( ioQQQ, " Illegal value for thermionic efficiency in %s\n",chFile );
2862 fprintf( ioQQQ, " Line #%ld, value=%.4e\n",dl,gd->therm_eff);
2863 fprintf( ioQQQ, " Allowed values are 0. < efficiency <= 1.\n");
2865 }
2866
2867 /* sublimation temperature in K */
2868 mie_next_data(chFile,io2,chLine,&dl);
2869 mie_read_double(chFile,chLine,&gd->subl_temp,true,dl);
2870
2871 /* >>chng 02 sep 18, add keyword for special files (grey grains, PAH's, etc.), PvH */
2872 mie_next_data(chFile,io2,chLine,&dl);
2873 mie_read_word(chLine,chWord,WORDLEN,true);
2874
2875 if( nMatch( "RFI_", chWord ) )
2876 gd->rfiType = RFI_TABLE;
2877 else if( nMatch( "OPC_", chWord ) )
2878 gd->rfiType = OPC_TABLE;
2879 else if( nMatch( "GREY", chWord ) )
2880 gd->rfiType = OPC_GREY;
2881 else if( nMatch( "PAH1", chWord ) )
2882 gd->rfiType = OPC_PAH1;
2883 else if( nMatch( "PH2N", chWord ) )
2884 gd->rfiType = OPC_PAH2N;
2885 else if( nMatch( "PH2C", chWord ) )
2886 gd->rfiType = OPC_PAH2C;
2887 else if( nMatch( "PH3N", chWord ) )
2888 gd->rfiType = OPC_PAH3N;
2889 else if( nMatch( "PH3C", chWord ) )
2890 gd->rfiType = OPC_PAH3C;
2891 else
2892 {
2893 fprintf( ioQQQ, " Illegal keyword in %s\n",chFile );
2894 fprintf( ioQQQ, " Line #%ld, value=%s\n",dl,chWord);
2895 fprintf( ioQQQ, " Allowed values are: RFI_TBL, OPC_TBL, GREY, PAH1, PH2N, PH2C, PH3N, PH3C\n");
2897 }
2898
2899 switch( gd->rfiType )
2900 {
2901 case RFI_TABLE:
2902 /* nridf is for choosing ref index or diel function input
2903 * case 2 allows greater accuracy reading in, when nr is close to 1. */
2904 mie_next_data(chFile,io2,chLine,&dl);
2905 mie_read_long(chFile,chLine,&nridf,true,dl);
2906 if( nridf > 3 )
2907 {
2908 fprintf( ioQQQ, " Illegal data code in %s\n",chFile );
2909 fprintf( ioQQQ, " Line #%ld, data code=%ld\n",dl,nridf);
2911 }
2912
2913 /* no. of principal axes, always 1 for amorphous grains,
2914 * maybe larger for crystalline grains */
2915 mie_next_data(chFile,io2,chLine,&dl);
2916 mie_read_long(chFile,chLine,&gd->nAxes,true,dl);
2917 if( gd->nAxes > NAX )
2918 {
2919 fprintf( ioQQQ, " Illegal no. of axes in %s\n",chFile );
2920 fprintf( ioQQQ, " Line #%ld, number=%ld\n",dl,gd->nAxes);
2922 }
2923
2924 /* now get relative weights of axes */
2925 mie_next_data(chFile,io2,chLine,&dl);
2926 switch( gd->nAxes )
2927 {
2928 case 1:
2929 mie_read_double(chFile,chLine,&gd->wt[0],true,dl);
2930 total = gd->wt[0];
2931 break;
2932 case 2:
2933 if( sscanf( chLine, "%lf %lf", &gd->wt[0], &gd->wt[1] ) != 2 )
2934 {
2935 fprintf( ioQQQ, " Syntax error in %s\n",chFile);
2936 fprintf( ioQQQ, " Line #%ld: %s\n",dl,chLine);
2938 }
2939 if( gd->wt[0] <= 0. || gd->wt[1] <= 0. )
2940 {
2941 fprintf( ioQQQ, " Illegal data in %s\n",chFile);
2942 fprintf( ioQQQ, " Line #%ld: %s\n",dl,chLine);
2944 }
2945 total = gd->wt[0] + gd->wt[1];
2946 break;
2947 case 3:
2948 if( sscanf( chLine, "%lf %lf %lf", &gd->wt[0], &gd->wt[1], &gd->wt[2] ) != 3 )
2949 {
2950 fprintf( ioQQQ, " Syntax error in %s\n",chFile);
2951 fprintf( ioQQQ, " Line #%ld: %s\n",dl,chLine);
2953 }
2954 if( gd->wt[0] <= 0. || gd->wt[1] <= 0. || gd->wt[2] <= 0. )
2955 {
2956 fprintf( ioQQQ, " Illegal data in %s\n",chFile);
2957 fprintf( ioQQQ, " Line #%ld: %s\n",dl,chLine);
2959 }
2960 total = gd->wt[0] + gd->wt[1] + gd->wt[2];
2961 break;
2962 default:
2963 fprintf( ioQQQ, " insane case for gd->nAxes: %ld\n", gd->nAxes );
2964 ShowMe();
2966 }
2967 for( j=0; j < gd->nAxes; j++ )
2968 {
2969 gd->wt[j] /= total;
2970
2971 /* read in optical constants for each principal axis. */
2972 mie_next_data(chFile,io2,chLine,&dl);
2973 mie_read_long(chFile,chLine,&gd->ndata[j],false,dl);
2974 if( gd->ndata[j] < 2 )
2975 {
2976 fprintf( ioQQQ, " Illegal number of data points in %s\n",chFile );
2977 fprintf( ioQQQ, " Line #%ld, number=%ld\n",dl,gd->ndata[j]);
2979 }
2980
2981 /* allocate space for wavelength and optical constants arrays */
2982 gd->wavlen[j].resize( gd->ndata[j] );
2983 gd->n[j].resize( gd->ndata[j] );
2984 gd->nr1[j].resize( gd->ndata[j] );
2985
2986 for( i=0; i < gd->ndata[j]; i++ )
2987 {
2988 /* read in the wavelength in microns
2989 * and the complex refractive index or dielectric function of material */
2990 mie_next_data(chFile,io2,chLine,&dl);
2991 if( sscanf( chLine, "%lf %lf %lf", &gd->wavlen[j][i], &nr, &ni ) != 3 )
2992 {
2993 fprintf( ioQQQ, " Syntax error in %s\n",chFile);
2994 fprintf( ioQQQ, " Line #%ld: %s\n",dl,chLine);
2996 }
2997 if( gd->wavlen[j][i] <= 0. )
2998 {
2999 fprintf( ioQQQ, " Illegal value for wavelength in %s\n",chFile);
3000 fprintf( ioQQQ, " Line #%ld, wavl=%14.6e\n",dl,gd->wavlen[j][i]);
3002 }
3003 /* the data in the input file should be sorted on wavelength, either
3004 * strictly monotonically increasing or decreasing, check this here... */
3005 if( i == 1 )
3006 {
3007 sgn = sign3(gd->wavlen[j][1]-gd->wavlen[j][0]);
3008 if( sgn == 0 )
3009 {
3010 fprintf( ioQQQ, " Illegal value for wavelength in %s\n",chFile);
3011 fprintf( ioQQQ, " Line #%ld, wavl=%14.6e\n",dl,gd->wavlen[j][i]);
3013 }
3014 }
3015 else if( i > 1 )
3016 {
3017 if( sign3(gd->wavlen[j][i]-gd->wavlen[j][i-1]) != sgn )
3018 {
3019 fprintf( ioQQQ, " Illegal value for wavelength in %s\n",chFile);
3020 fprintf( ioQQQ, " Line #%ld, wavl=%14.6e\n",dl,gd->wavlen[j][i]);
3022 }
3023 }
3024 /* this version reads in real and imaginary parts of the refractive
3025 * index, with imaginary part positive (nridf = 3) or nr-1 (nridf = 2) or
3026 * real and imaginary parts of the dielectric function (nridf = 1) */
3027 switch( nridf )
3028 {
3029 case 1:
3030 eps1 = nr;
3031 eps2 = ni;
3032 dftori(&nr,&ni,eps1,eps2);
3033 gd->nr1[j][i] = nr - 1.;
3034 break;
3035 case 2:
3036 gd->nr1[j][i] = nr;
3037 nr += 1.;
3038 break;
3039 case 3:
3040 gd->nr1[j][i] = nr - 1.;
3041 break;
3042 default:
3043 fprintf( ioQQQ, " insane case for nridf: %ld\n", nridf );
3044 ShowMe();
3046 }
3047 gd->n[j][i] = complex<double>(nr,ni);
3048
3049 /* sanity checks */
3050 if( nr <= 0. || ni < 0. )
3051 {
3052 fprintf( ioQQQ, " Illegal value for refractive index in %s\n",chFile);
3053 fprintf( ioQQQ, " Line #%ld, (nr,ni)=(%14.6e,%14.6e)\n",dl,nr,ni);
3055 }
3056 ASSERT( fabs(nr-1.-gd->nr1[j][i]) < 10.*nr*DBL_EPSILON );
3057 }
3058 }
3059 break;
3060 case OPC_TABLE:
3061 /* no. of data columns in OPC_TABLE file:
3062 * 1: absorption cross sections only
3063 * 2: absorption + scattering cross sections
3064 * 3: absorption + pure scattering cross sections + asymmetry factor
3065 * 4: absorption + pure scattering cross sections + asymmetry factor +
3066 * inverse attenuation length */
3067 mie_next_data(chFile,io2,chLine,&dl);
3068 mie_read_long(chFile,chLine,&gd->nOpcCols,true,dl);
3069 if( gd->nOpcCols > NDAT )
3070 {
3071 fprintf( ioQQQ, " Illegal no. of data columns in %s\n",chFile );
3072 fprintf( ioQQQ, " Line #%ld, number=%ld\n",dl,gd->nOpcCols);
3074 }
3075
3076 /* keyword indicating whether the table contains linear or logarithmic data */
3077 mie_next_data(chFile,io2,chLine,&dl);
3078 mie_read_word(chLine,chWord,WORDLEN,true);
3079
3080 if( nMatch( "LINE", chWord ) )
3081 lgLogData = false;
3082 else if( nMatch( "LOG", chWord ) )
3083 lgLogData = true;
3084 else
3085 {
3086 fprintf( ioQQQ, " Keyword not recognized in %s\n",chFile );
3087 fprintf( ioQQQ, " Line #%ld, keyword=%s\n",dl,chWord);
3089 }
3090
3091
3092 /* read in number of data points supplied. */
3093 mie_next_data(chFile,io2,chLine,&dl);
3094 mie_read_long(chFile,chLine,&gd->nOpcData,false,dl);
3095 if( gd->nOpcData < 2 )
3096 {
3097 fprintf( ioQQQ, " Illegal number of data points in %s\n",chFile );
3098 fprintf( ioQQQ, " Line #%ld, number=%ld\n",dl,gd->nOpcData);
3100 }
3101
3102 /* allocate space for frequency and data arrays */
3103 gd->opcAnu.resize(gd->nOpcData);
3104 for( j=0; j < gd->nOpcCols; j++ )
3105 gd->opcData[j].resize(gd->nOpcData);
3106
3107 tmp1 = -log10(1.01*DBL_MIN);
3108 tmp2 = log10(0.99*DBL_MAX);
3109 LargestLog = min(tmp1,tmp2);
3110
3111 /* now read the data... each line should contain:
3112 *
3113 * if gd->nOpcCols == 1, anu, abs_cs
3114 * if gd->nOpcCols == 2, anu, abs_cs, sct_cs
3115 * if gd->nOpcCols == 3, anu, abs_cs, sct_cs, inv_att_len
3116 *
3117 * the frequencies in the table should be either monotonically increasing or decreasing.
3118 * frequencies should be in Ryd, cross sections in cm^2/H, and the inverse attenuation length
3119 * in cm^-1. If lgLogData is true, each number should be the log10 of the data value */
3120 for( i=0; i < gd->nOpcData; i++ )
3121 {
3122 mie_next_data(chFile,io2,chLine,&dl);
3123 switch( gd->nOpcCols )
3124 {
3125 case 1:
3126 if( sscanf( chLine, "%lf %lf", &gd->opcAnu[i], &gd->opcData[0][i] ) != 2 )
3127 {
3128 fprintf( ioQQQ, " Syntax error in %s\n",chFile);
3129 fprintf( ioQQQ, " Line #%ld: %s\n",dl,chLine);
3131 }
3132 break;
3133 case 2:
3134 if( sscanf( chLine, "%lf %lf %lf", &gd->opcAnu[i], &gd->opcData[0][i],
3135 &gd->opcData[1][i] ) != 3 )
3136 {
3137 fprintf( ioQQQ, " Syntax error in %s\n",chFile);
3138 fprintf( ioQQQ, " Line #%ld: %s\n",dl,chLine);
3140 }
3141 break;
3142 case 3:
3143 if( sscanf( chLine, "%lf %lf %lf %lf", &gd->opcAnu[i], &gd->opcData[0][i],
3144 &gd->opcData[1][i], &gd->opcData[2][i] ) != 4 )
3145 {
3146 fprintf( ioQQQ, " Syntax error in %s\n",chFile);
3147 fprintf( ioQQQ, " Line #%ld: %s\n",dl,chLine);
3149 }
3150 break;
3151 case 4:
3152 if( sscanf( chLine, "%lf %lf %lf %lf %lf", &gd->opcAnu[i], &gd->opcData[0][i],
3153 &gd->opcData[1][i], &gd->opcData[2][i], &gd->opcData[3][i] ) != 5 )
3154 {
3155 fprintf( ioQQQ, " Syntax error in %s\n",chFile);
3156 fprintf( ioQQQ, " Line #%ld: %s\n",dl,chLine);
3158 }
3159 break;
3160 default:
3161 fprintf( ioQQQ, " insane case for gd->nOpcCols: %ld\n", gd->nOpcCols );
3162 ShowMe();
3164 }
3165 if( lgLogData )
3166 {
3167 /* this test will guarantee there will be neither under- nor overflows */
3168 if( fabs(gd->opcAnu[i]) > LargestLog )
3169 {
3170 fprintf( ioQQQ, " Illegal value for frequency in %s\n",chFile );
3171 fprintf( ioQQQ, " Line #%ld, freq=%14.6e\n",dl,gd->opcAnu[i] );
3173 }
3174 gd->opcAnu[i] = pow(10.,gd->opcAnu[i]);
3175 for( j=0; j < gd->nOpcCols; j++ )
3176 {
3177 if( fabs(gd->opcData[j][i]) > LargestLog )
3178 {
3179 fprintf( ioQQQ, " Illegal data value in %s\n",chFile );
3180 fprintf( ioQQQ, " Line #%ld, value=%14.6e\n",dl,gd->opcData[j][i] );
3182 }
3183 gd->opcData[j][i] = pow(10.,gd->opcData[j][i]);
3184 }
3185 }
3186 if( gd->opcAnu[i] <= 0. )
3187 {
3188 fprintf( ioQQQ, " Illegal value for frequency in %s\n",chFile );
3189 fprintf( ioQQQ, " Line #%ld, freq=%14.6e\n",dl,gd->opcAnu[i] );
3191 }
3192 for( j=0; j < gd->nOpcCols; j++ )
3193 {
3194 if( gd->opcData[j][i] <= 0. )
3195 {
3196 fprintf( ioQQQ, " Illegal data value in %s\n",chFile );
3197 fprintf( ioQQQ, " Line #%ld, value=%14.6e\n",dl,gd->opcData[j][i] );
3199 }
3200 }
3201 /* the data in the input file should be sorted on frequency, either
3202 * strictly monotonically increasing or decreasing, check this here... */
3203 if( i == 1 )
3204 {
3205 sgn = sign3(gd->opcAnu[1]-gd->opcAnu[0]);
3206 if( sgn == 0 )
3207 {
3208 double dataVal = lgLogData ? log10(gd->opcAnu[1]) : gd->opcAnu[1];
3209 fprintf( ioQQQ, " Illegal value for frequency in %s\n",chFile );
3210 fprintf( ioQQQ, " Line #%ld, freq=%14.6e\n",dl,dataVal );
3212 }
3213 }
3214 else if( i > 1 )
3215 {
3216 if( sign3(gd->opcAnu[i]-gd->opcAnu[i-1]) != sgn )
3217 {
3218 double dataVal = lgLogData ? log10(gd->opcAnu[i]) : gd->opcAnu[i];
3219 fprintf( ioQQQ, " Illegal value for frequency in %s\n",chFile);
3220 fprintf( ioQQQ, " Line #%ld, freq=%14.6e\n",dl,dataVal);
3222 }
3223 }
3224 }
3225 gd->nAxes = 1;
3226 break;
3227 case OPC_GREY:
3228 case OPC_PAH1:
3229 case OPC_PAH2N:
3230 case OPC_PAH2C:
3231 case OPC_PAH3N:
3232 case OPC_PAH3C:
3233 /* nothing much to be done here, the opacities
3234 * will be calculated without any further data */
3235 gd->nAxes = 1;
3236 break;
3237 default:
3238 fprintf( ioQQQ, " Insane value for gd->rfiType: %d, bailing out\n", gd->rfiType );
3239 ShowMe();
3241 }
3242
3243 fclose(io2);
3244 return;
3245}
3246
3247
3248/* construct optical constants for mixed grain using a specific EMT */
3249STATIC void mie_read_mix(/*@in@*/ const char *chFile,
3250 /*@out@*/ grain_data *gd)
3251{
3252 emt_type EMTtype;
3253 long int dl,
3254 i,
3255 j,
3256 k,
3257 l,
3258 nelem,
3259 nMaterial,
3260 sumAxes;
3261 double maxIndex = DBL_MAX,
3262 minIndex = DBL_MAX,
3263 nAtoms,
3264 sum,
3265 sum2,
3266 wavHi,
3267 wavLo,
3268 wavStep;
3269 complex<double> eps_eff(-DBL_MAX,-DBL_MAX);
3270 char chLine[FILENAME_PATH_LENGTH_2],
3271 chWord[FILENAME_PATH_LENGTH_2],
3272 *str;
3273 FILE *io2;
3274
3275 DEBUG_ENTRY( "mie_read_mix()" );
3276
3277 io2 = open_data( chFile, "r", AS_LOCAL_ONLY );
3278
3279 dl = 0; /* line counter for input file */
3280
3281 /* first read magic number */
3282 mie_next_data(chFile,io2,chLine,&dl);
3283 mie_read_long(chFile,chLine,&gd->magic,true,dl);
3284 if( gd->magic != MAGIC_MIX )
3285 {
3286 fprintf( ioQQQ, " Mixed grain file %s has obsolete magic number\n",chFile );
3287 fprintf( ioQQQ, " I found %ld, but expected %ld on line #%ld\n",gd->magic,MAGIC_MIX,dl );
3288 fprintf( ioQQQ, " Please replace this file with an up to date version\n" );
3290 }
3291
3292 /* default depletion of grain molecule */
3293 mie_next_data(chFile,io2,chLine,&dl);
3294 mie_read_double(chFile,chLine,&gd->depl,true,dl);
3295 if( gd->depl > 1. )
3296 {
3297 fprintf( ioQQQ, " Illegal value for default depletion in %s\n",chFile );
3298 fprintf( ioQQQ, " Line #%ld, depl=%14.6e\n",dl,gd->depl);
3300 }
3301
3302 /* read number of different materials contained in this grain */
3303 mie_next_data(chFile,io2,chLine,&dl);
3304 mie_read_long(chFile,chLine,&nMaterial,true,dl);
3305 if( nMaterial < 2 )
3306 {
3307 fprintf( ioQQQ, " Illegal number of materials in mixed grain file %s\n",chFile );
3308 fprintf( ioQQQ, " I found %ld on line #%ld\n",nMaterial,dl );
3309 fprintf( ioQQQ, " This number should be at least 2\n" );
3311 }
3312
3313 vector<double> frac(nMaterial);
3314 vector<double> frac2(nMaterial);
3315 vector<grain_data> gdArr(nMaterial);
3316
3317 sum = 0.;
3318 sum2 = 0.;
3319 sumAxes = 0;
3320 for( i=0; i < nMaterial; i++ )
3321 {
3322 char chFile2[FILENAME_PATH_LENGTH_2];
3323
3324 /* each line contains relative fraction of volume occupied by each material,
3325 * followed by the name of the refractive index file between double quotes */
3326 mie_next_data(chFile,io2,chLine,&dl);
3327 mie_read_double(chFile,chLine,&frac[i],true,dl);
3328
3329 sum += frac[i];
3330
3331 str = strchr_s( chLine, '\"' );
3332 if( str != NULL )
3333 {
3334 strcpy( chFile2, ++str );
3335 str = strchr_s( chFile2, '\"' );
3336 if( str != NULL )
3337 *str = '\0';
3338 }
3339 if( str == NULL )
3340 {
3341 fprintf( ioQQQ, " No pair of double quotes was found on line #%ld of file %s\n",dl,chFile );
3342 fprintf( ioQQQ, " Please supply the refractive index file name between double quotes\n" );
3344 }
3345
3346 mie_read_rfi( chFile2, &gdArr[i] );
3347 if( gdArr[i].rfiType != RFI_TABLE )
3348 {
3349 fprintf( ioQQQ, " Input error on line #%ld of file %s\n",dl,chFile );
3350 fprintf( ioQQQ, " File %s is not of type RFI_TBL, this is illegal\n",chFile2 );
3352 }
3353
3354 /* this test also guarantees that sum2 > 0 */
3355 if( i == nMaterial-1 && gdArr[i].mol_weight == 0. )
3356 {
3357 fprintf( ioQQQ, " Input error on line #%ld of file %s\n",dl,chFile );
3358 fprintf( ioQQQ, " Last entry in list of materials is vacuum, this is not allowed\n" );
3359 fprintf( ioQQQ, " Please move this entry to an earlier position in the file\n" );
3361 }
3362
3363 frac2[i] = ( gdArr[i].mol_weight > 0. ) ? frac[i] : 0.;
3364 sum2 += frac2[i];
3365
3366 sumAxes += gdArr[i].nAxes;
3367 }
3368
3369 /* read keyword to determine which EMT to use */
3370 mie_next_data(chFile,io2,chLine,&dl);
3371 mie_read_word(chLine,chWord,WORDLEN,true);
3372
3373 if( nMatch( "FA00", chWord ) )
3374 EMTtype = FARAFONOV00;
3375 else if( nMatch( "ST95", chWord ) )
3376 EMTtype = STOGNIENKO95;
3377 else if( nMatch( "BR35", chWord ) )
3378 EMTtype = BRUGGEMAN35;
3379 else
3380 {
3381 fprintf( ioQQQ, " Keyword not recognized in %s\n",chFile );
3382 fprintf( ioQQQ, " Line #%ld, keyword=%s\n",dl,chWord);
3384 }
3385
3386 /* normalize sum of fractional volumes to 1 */
3387 for( i=0; i < nMaterial; i++ )
3388 {
3389 frac[i] /= sum;
3390 frac2[i] /= sum2;
3391 /* renormalize elmAbun to chemical formula */
3392 for( nelem=0; nelem < LIMELM; nelem++ )
3393 {
3394 gdArr[i].elmAbun[nelem] /= gdArr[i].abun*gdArr[i].depl;
3395 }
3396 }
3397
3398 wavLo = 0.;
3399 wavHi = DBL_MAX;
3400 gd->abun = DBL_MAX;
3401 for( nelem=0; nelem < LIMELM; nelem++ )
3402 {
3403 gd->elmAbun[nelem] = 0.;
3404 }
3405 gd->mol_weight = 0.;
3406 gd->rho = 0.;
3407 gd->work = DBL_MAX;
3408 gd->bandgap = DBL_MAX;
3409 gd->therm_eff = 0.;
3410 gd->subl_temp = DBL_MAX;
3411 gd->nAxes = 1;
3412 gd->wt[0] = 1.;
3413 gd->ndata[0] = MIX_TABLE_SIZE;
3414 gd->rfiType = RFI_TABLE;
3415
3416 for( i=0; i < nMaterial; i++ )
3417 {
3418 for( k=0; k < gdArr[i].nAxes; k++ )
3419 {
3420 double wavMin = min(gdArr[i].wavlen[k][0],gdArr[i].wavlen[k][gdArr[i].ndata[k]-1]);
3421 double wavMax = max(gdArr[i].wavlen[k][0],gdArr[i].wavlen[k][gdArr[i].ndata[k]-1]);
3422 wavLo = max(wavLo,wavMin);
3423 wavHi = min(wavHi,wavMax);
3424 }
3425 minIndex = DBL_MAX;
3426 maxIndex = 0.;
3427 for( nelem=0; nelem < LIMELM; nelem++ )
3428 {
3429 gd->elmAbun[nelem] += frac[i]*gdArr[i].elmAbun[nelem];
3430 if( gd->elmAbun[nelem] > 0. )
3431 {
3432 minIndex = min(minIndex,gd->elmAbun[nelem]);
3433 }
3434 maxIndex = max(maxIndex,gd->elmAbun[nelem]);
3435 }
3436 gd->mol_weight += frac[i]*gdArr[i].mol_weight;
3437 gd->rho += frac[i]*gdArr[i].rho;
3438 /* ignore parameters for vacuum */
3439 if( gdArr[i].mol_weight > 0. )
3440 {
3441 gd->abun = min(gd->abun,gdArr[i].abun/frac[i]);
3442 switch( EMTtype )
3443 {
3444 case FARAFONOV00:
3445 /* this is appropriate for a layered grain */
3446 gd->work = gdArr[i].work;
3447 gd->bandgap = gdArr[i].bandgap;
3448 gd->therm_eff = gdArr[i].therm_eff;
3449 break;
3450 case STOGNIENKO95:
3451 case BRUGGEMAN35:
3452 /* this is appropriate for a randomly mixed grain */
3453 gd->work = min(gd->work,gdArr[i].work);
3454 gd->bandgap = min(gd->bandgap,gdArr[i].bandgap);
3455 gd->therm_eff += frac2[i]*gdArr[i].therm_eff;
3456 break;
3457 default:
3458 fprintf( ioQQQ, " Insanity in mie_read_mix\n" );
3459 ShowMe();
3461 }
3462 gd->matType = gdArr[i].matType;
3463 gd->subl_temp = min(gd->subl_temp,gdArr[i].subl_temp);
3464 }
3465 }
3466
3467 if( gd->rho <= 0. )
3468 {
3469 fprintf( ioQQQ, " Illegal value for the density: %.3e\n", gd->rho );
3471 }
3472 if( gd->mol_weight <= 0. )
3473 {
3474 fprintf( ioQQQ, " Illegal value for the molecular weight: %.3e\n", gd->mol_weight );
3476 }
3477 if( maxIndex <= 0. )
3478 {
3479 fprintf( ioQQQ, " No atoms were found in the grain molecule\n" );
3481 }
3482
3483 /* further sanity checks */
3484 ASSERT( wavLo > 0. && wavHi < DBL_MAX && wavLo < wavHi );
3485 ASSERT( gd->abun > 0. && gd->abun < DBL_MAX );
3486 ASSERT( gd->work > 0. && gd->work < DBL_MAX );
3487 ASSERT( gd->bandgap >= 0. && gd->bandgap < gd->work );
3488 ASSERT( gd->therm_eff > 0. && gd->therm_eff <= 1. );
3489 ASSERT( gd->subl_temp > 0. && gd->subl_temp < DBL_MAX );
3490 ASSERT( minIndex > 0. && minIndex < DBL_MAX );
3491
3492 /* apply safety margin */
3493 wavLo *= 1. + 10.*DBL_EPSILON;
3494 wavHi *= 1. - 10.*DBL_EPSILON;
3495
3496 /* renormalize the chemical formula such that the lowest index is 1 */
3497 nAtoms = 0.;
3498 for( nelem=0; nelem < LIMELM; nelem++ )
3499 {
3500 gd->elmAbun[nelem] /= minIndex;
3501 nAtoms += gd->elmAbun[nelem];
3502 }
3503 ASSERT( nAtoms > 0. );
3504 gd->abun *= minIndex;
3505 gd->mol_weight /= minIndex;
3506 /* calculate average weight per atom */
3507 gd->atom_weight = gd->mol_weight/nAtoms;
3508
3509 mie_write_form(gd->elmAbun,chWord);
3510 fprintf( ioQQQ, "\n The chemical formula of the new grain molecule is: %s\n", chWord );
3511 fprintf( ioQQQ, " The abundance wrt H at maximum depletion of this molecule is: %.3e\n",
3512 gd->abun );
3513 fprintf( ioQQQ, " The abundance wrt H at standard depletion of this molecule is: %.3e\n",
3514 gd->abun*gd->depl );
3515
3516 /* finally renormalize elmAbun back to abundance relative to hydrogen */
3517 for( nelem=0; nelem < LIMELM; nelem++ )
3518 {
3519 gd->elmAbun[nelem] *= gd->abun*gd->depl;
3520 }
3521
3522 vector<double> delta(sumAxes);
3523 vector<double> frdelta(sumAxes);
3524 vector< complex<double> > eps(sumAxes);
3525
3526 l = 0;
3527 for( i=0; i < nMaterial; i++ )
3528 {
3529 for( k=0; k < gdArr[i].nAxes; k++ )
3530 {
3531 frdelta[l] = gdArr[i].wt[k]*frac[i];
3532 delta[l] = ( l == 0 ) ? frdelta[l] : delta[l-1] + frdelta[l];
3533 ++l;
3534 }
3535 }
3536 ASSERT( l == sumAxes && fabs(delta[l-1]-1.) < 10.*DBL_EPSILON );
3537
3538 /* allocate space for wavelength and optical constants arrays */
3539 gd->wavlen[0].resize( gd->ndata[0] );
3540 gd->n[0].resize( gd->ndata[0] );
3541 gd->nr1[0].resize( gd->ndata[0] );
3542
3543 wavStep = log(wavHi/wavLo)/(double)(gd->ndata[0]-1);
3544
3545 switch( EMTtype )
3546 {
3547 case FARAFONOV00:
3548 /* this implements the EMT described in
3549 * >>refer grain physics Voshchinnikov N.V., Mathis J.S., 1999, ApJ, 526, 257
3550 * based on the theory described in
3551 * >>refer grain physics Farafonov V.G., 2000, Optics & Spectroscopy, 88, 441
3552 *
3553 * NB - note that Eq. 3 in Voshchinnikov & Mathis is incorrect! */
3554 for( j=0; j < gd->ndata[0]; j++ )
3555 {
3556 double nre,nim;
3557 complex<double> a1,a2,a1c,a2c,a11,a12,a21,a22,ratio;
3558
3559 gd->wavlen[0][j] = wavLo*exp((double)j*wavStep);
3560
3561 init_eps(gd->wavlen[0][j],nMaterial,gdArr,eps);
3562
3563 ratio = eps[0]/eps[1];
3564
3565 a1 = (ratio+2.)/3.;
3566 a2 = (1.-ratio)*delta[0];
3567
3568 for( l=1; l < sumAxes-1; l++ )
3569 {
3570 ratio = eps[l]/eps[l+1];
3571
3572 a1c = a1;
3573 a2c = a2;
3574 a11 = (ratio+2.)/3.;
3575 a12 = (2.-2.*ratio)/(9.*delta[l]);
3576 a21 = (1.-ratio)*delta[l];
3577 a22 = (2.*ratio+1.)/3.;
3578
3579 a1 = a11*a1c + a12*a2c;
3580 a2 = a21*a1c + a22*a2c;
3581 }
3582
3583 a1c = a1;
3584 a2c = a2;
3585 a11 = 1.;
3586 a12 = 1./3.;
3587 a21 = eps[sumAxes-1];
3588 a22 = -2./3.*eps[sumAxes-1];
3589
3590 a1 = a11*a1c + a12*a2c;
3591 a2 = a21*a1c + a22*a2c;
3592
3593 ratio = a2/a1;
3594 dftori(&nre,&nim,ratio.real(),ratio.imag());
3595
3596 gd->n[0][j] = complex<double>(nre,nim);
3597 gd->nr1[0][j] = nre-1.;
3598 }
3599 break;
3600 case STOGNIENKO95:
3601 case BRUGGEMAN35:
3602 for( j=0; j < gd->ndata[0]; j++ )
3603 {
3604 const double EPS_TOLER = 10.*DBL_EPSILON;
3605 double nre,nim;
3606 complex<double> eps0;
3607
3608 gd->wavlen[0][j] = wavLo*exp((double)j*wavStep);
3609
3610 init_eps(gd->wavlen[0][j],nMaterial,gdArr,eps);
3611
3612 /* get initial estimate for effective dielectric function */
3613 if( j == 0 )
3614 {
3615 /* use simple average as first estimate */
3616 eps0 = 0.;
3617 for( l=0; l < sumAxes; l++ )
3618 eps0 += frdelta[l]*eps[l];
3619 }
3620 else
3621 {
3622 /* use solution from previous wavelength as first estimate */
3623 eps0 = eps_eff;
3624 }
3625
3626 if( EMTtype == STOGNIENKO95 )
3627 /* this implements the EMT described in
3628 * >>refer grain physics Stognienko R., Henning Th., Ossenkopf V., 1995, A&A, 296, 797 */
3629 eps_eff = cnewton( Stognienko, frdelta, eps, sumAxes, eps0, EPS_TOLER );
3630 else if( EMTtype == BRUGGEMAN35 )
3631 /* this implements the classical Bruggeman rule
3632 * >>refer grain physics Bruggeman D.A.G., 1935, Ann. Phys. (5th series), 24, 636 */
3633 eps_eff = cnewton( Bruggeman, frdelta, eps, sumAxes, eps0, EPS_TOLER );
3634 else
3635 {
3636 fprintf( ioQQQ, " Insanity in mie_read_mix\n" );
3637 ShowMe();
3639 }
3640
3641 dftori(&nre,&nim,eps_eff.real(),eps_eff.imag());
3642
3643 gd->n[0][j] = complex<double>(nre,nim);
3644 gd->nr1[0][j] = nre-1.;
3645 }
3646 break;
3647 default:
3648 fprintf( ioQQQ, " Insanity in mie_read_mix\n" );
3649 ShowMe();
3651 }
3652 return;
3653}
3654
3655
3656/* helper routine for mie_read_mix, initializes the array of dielectric functions */
3657STATIC void init_eps(double wavlen,
3658 long nMaterial,
3659 /*@in@*/ const vector<grain_data>& gdArr, /* gdArr[nMaterial] */
3660 /*@out@*/ vector< complex<double> >& eps) /* eps[sumAxes] */
3661{
3662 long i,
3663 k;
3664
3665 long l = 0;
3666
3667 DEBUG_ENTRY( "init_eps()" );
3668 for( i=0; i < nMaterial; i++ )
3669 {
3670 for( k=0; k < gdArr[i].nAxes; k++ )
3671 {
3672 bool lgErr;
3673 long ind;
3674 double eps1,eps2,frc,nim,nre;
3675
3676 find_arr(wavlen,gdArr[i].wavlen[k],gdArr[i].ndata[k],&ind,&lgErr);
3677 ASSERT( !lgErr );
3678 frc = (wavlen-gdArr[i].wavlen[k][ind])/(gdArr[i].wavlen[k][ind+1]-gdArr[i].wavlen[k][ind]);
3679 ASSERT( frc > 0.-10.*DBL_EPSILON && frc < 1.+10.*DBL_EPSILON );
3680 nre = (1.-frc)*gdArr[i].n[k][ind].real() + frc*gdArr[i].n[k][ind+1].real();
3681 ASSERT( nre > 0. );
3682 nim = (1.-frc)*gdArr[i].n[k][ind].imag() + frc*gdArr[i].n[k][ind+1].imag();
3683 ASSERT( nim >= 0. );
3684 ritodf(nre,nim,&eps1,&eps2);
3685 eps[l++] = complex<double>(eps1,eps2);
3686 }
3687 }
3688 return;
3689}
3690
3691
3692/*******************************************************
3693 * This routine is derived from the routine Znewton *
3694 * --------------------------------------------------- *
3695 * Reference; BASIC Scientific Subroutines, Vol. II *
3696 * by F.R. Ruckdeschel, BYTE/McGRAWW-HILL, 1981. *
3697 * *
3698 * C++ version by J-P Moreau, Paris. *
3699 ******************************************************/
3700/* find complex root of fun using the Newton-Raphson algorithm */
3701STATIC complex<double> cnewton(
3702 void(*fun)(complex<double>,const vector<double>&,const vector< complex<double> >&,
3703 long,complex<double>*,double*,double*),
3704 /*@in@*/ const vector<double>& frdelta, /* frdelta[sumAxes] */
3705 /*@in@*/ const vector< complex<double> >& eps, /* eps[sumAxes] */
3706 long sumAxes,
3707 complex<double> x0,
3708 double tol)
3709{
3710 int i;
3711
3712 const int LOOP_MAX = 100;
3713 const double TINY = 1.e-12;
3714
3715 DEBUG_ENTRY( "cnewton()" );
3716 for( i=0; i < LOOP_MAX; i++ )
3717 {
3718 complex<double> x1,y;
3719 double dudx,dudy,norm2;
3720
3721 (*fun)(x0,frdelta,eps,sumAxes,&y,&dudx,&dudy);
3722
3723 norm2 = pow2(dudx) + pow2(dudy);
3724 /* guard against norm2 == 0 */
3725 if( norm2 < TINY*norm(y) )
3726 {
3727 fprintf( ioQQQ, " cnewton - zero divide error\n" );
3728 ShowMe();
3730 }
3731 x1 = x0 - complex<double>( y.real()*dudx-y.imag()*dudy, y.imag()*dudx+y.real()*dudy )/norm2;
3732
3733 /* check for convergence */
3734 if( fabs(x0.real()/x1.real()-1.) + fabs(x0.imag()/x1.imag()-1.) < tol )
3735 {
3736 return x1;
3737 }
3738
3739 x0 = x1;
3740 }
3741
3742 fprintf( ioQQQ, " cnewton did not converge\n" );
3743 ShowMe();
3745}
3746
3747
3748/* this evaluates the function defined in Eq. 3 of
3749 * >>refer grain physics Stognienko R., Henning Th., Ossenkopf V., 1995, A&A, 296, 797
3750 * and its derivatives */
3751STATIC void Stognienko(complex<double> x,
3752 /*@in@*/ const vector<double>& frdelta, /* frdelta[sumAxes] */
3753 /*@in@*/ const vector< complex<double> >& eps, /* eps[sumAxes] */
3754 long sumAxes,
3755 /*@out@*/ complex<double> *f,
3756 /*@out@*/ double *dudx,
3757 /*@out@*/ double *dudy)
3758{
3759 long i,
3760 l;
3761
3762 static const double L[4] = { 0., 1./2., 1., 1./3. };
3763 static const double fl[4] = { 5./9., 2./9., 2./9., 1. };
3764
3765 DEBUG_ENTRY( "Stognienko()" );
3766 *f = complex<double>(0.,0.);
3767 *dudx = 0.;
3768 *dudy = 0.;
3769 for( l=0; l < sumAxes; l++ )
3770 {
3771 complex<double> hlp = eps[l] - x;
3772 double h1 = eps[l].imag()*x.real() - eps[l].real()*x.imag();
3773
3774 for( i=0; i < 4; i++ )
3775 {
3776 double f1 = fl[i]*frdelta[l];
3777 double xx = ( i < 3 ) ? sin(PI*frdelta[l]) : cos(PI*frdelta[l]);
3778 complex<double> f2 = f1*xx*xx;
3779 complex<double> one = x + hlp*L[i];
3780 complex<double> two = f2*hlp/one;
3781 double h2 = norm(one);
3782 *f += two;
3783 *dudx -= f2.real()*(eps[l].real()*h2 + h1*2.*one.imag()*(1.-L[i]))/pow2(h2);
3784 *dudy -= f2.real()*(eps[l].imag()*h2 - h1*2.*one.real()*(1.-L[i]))/pow2(h2);
3785 }
3786 }
3787 return;
3788}
3789
3790
3791/* this evaluates the classical Bruggeman rule and its derivatives
3792 * >>refer grain physics Bruggeman D.A.G., 1935, Ann. Phys. (5th series), 24, 636 */
3793STATIC void Bruggeman(complex<double> x,
3794 /*@in@*/ const vector<double>& frdelta, /* frdelta[sumAxes] */
3795 /*@in@*/ const vector< complex<double> >& eps, /* eps[sumAxes] */
3796 long sumAxes,
3797 /*@out@*/ complex<double> *f,
3798 /*@out@*/ double *dudx,
3799 /*@out@*/ double *dudy)
3800{
3801 long l;
3802
3803 static const double L = 1./3.;
3804
3805 DEBUG_ENTRY( "Bruggeman()" );
3806 *f = complex<double>(0.,0.);
3807 *dudx = 0.;
3808 *dudy = 0.;
3809 for( l=0; l < sumAxes; l++ )
3810 {
3811 complex<double> hlp = eps[l] - x;
3812 double h1 = eps[l].imag()*x.real() - eps[l].real()*x.imag();
3813 complex<double> f2 = frdelta[l];
3814 complex<double> one = x + hlp*L;
3815 complex<double> two = f2*hlp/one;
3816 double h2 = norm(one);
3817 *f += two;
3818 *dudx -= f2.real()*(eps[l].real()*h2 + h1*2.*one.imag()*(1.-L))/pow2(h2);
3819 *dudy -= f2.real()*(eps[l].imag()*h2 - h1*2.*one.real()*(1.-L))/pow2(h2);
3820 }
3821 return;
3822}
3823
3824
3825/* read in the file with optical constants and other relevant information */
3826STATIC void mie_read_szd(/*@in@*/ const char *chFile,
3827 /*@out@*/ sd_data *sd)
3828{
3829 bool lgTryOverride = false;
3830 long int dl,
3831 i;
3832 double mul = 0.,
3833 ref_neg = DBL_MAX,
3834 ref_pos = DBL_MAX,
3835 step_neg = DBL_MAX,
3836 step_pos = DBL_MAX;
3837 char chLine[FILENAME_PATH_LENGTH_2],
3838 chWord[WORDLEN];
3839 FILE *io2;
3840
3841 /* these constants are used to get initial estimates for the cutoffs (lim)
3842 * in the SD_EXP_CUTOFFx and SD_xxx_NORMAL cases, they are iterated by
3843 * search_limit such that
3844 * lim^4 * dN/da(lim) == FRAC_CUTOFF * ref^4 * dN/da(ref)
3845 * where ref as an appropriate reference point for each of the cases */
3846 const double FRAC_CUTOFF = 1.e-4;
3847 const double MUL_CO1 = -log(FRAC_CUTOFF);
3848 const double MUL_CO2 = sqrt(MUL_CO1);
3849 const double MUL_CO3 = pow(MUL_CO1,1./3.);
3850 const double MUL_LND = sqrt(-2.*log(FRAC_CUTOFF));
3851 const double MUL_NRM = MUL_LND;
3852
3853 DEBUG_ENTRY( "mie_read_szd()" );
3854
3855 io2 = open_data( chFile, "r", AS_LOCAL_ONLY );
3856
3857 dl = 0; /* line counter for input file */
3858
3859 /* first read magic number */
3860 mie_next_data(chFile,io2,chLine,&dl);
3861 mie_read_long(chFile,chLine,&sd->magic,true,dl);
3862 if( sd->magic != MAGIC_SZD )
3863 {
3864 fprintf( ioQQQ, " Size distribution file %s has obsolete magic number\n",chFile );
3865 fprintf( ioQQQ, " I found %ld, but expected %ld on line #%ld\n",sd->magic,MAGIC_SZD,dl );
3866 fprintf( ioQQQ, " Please replace this file with an up to date version\n" );
3868 }
3869
3870 /* size distribution case */
3871 mie_next_data(chFile,io2,chLine,&dl);
3872 mie_read_word(chLine,chWord,WORDLEN,true);
3873
3874 if( nMatch( "SSIZ", chWord ) )
3875 {
3876 sd->sdCase = SD_SINGLE_SIZE;
3877 }
3878 else if( nMatch( "NCAR", chWord ) )
3879 {
3880 sd->sdCase = SD_NR_CARBON;
3881 }
3882 else if( nMatch( "POWE", chWord ) )
3883 {
3884 sd->sdCase = SD_POWERLAW;
3885 }
3886 else if( nMatch( "EXP1", chWord ) )
3887 {
3888 sd->sdCase = SD_EXP_CUTOFF1;
3889 sd->a[ipAlpha] = 1.;
3890 mul = MUL_CO1;
3891 }
3892 else if( nMatch( "EXP2", chWord ) )
3893 {
3894 sd->sdCase = SD_EXP_CUTOFF2;
3895 sd->a[ipAlpha] = 2.;
3896 mul = MUL_CO2;
3897 }
3898 else if( nMatch( "EXP3", chWord ) )
3899 {
3900 sd->sdCase = SD_EXP_CUTOFF3;
3901 sd->a[ipAlpha] = 3.;
3902 mul = MUL_CO3;
3903 }
3904 else if( nMatch( "LOGN", chWord ) )
3905 {
3906 sd->sdCase = SD_LOG_NORMAL;
3907 mul = MUL_LND;
3908 }
3909 /* this one must come after LOGNORMAL */
3910 else if( nMatch( "NORM", chWord ) )
3911 {
3912 sd->sdCase = SD_LIN_NORMAL;
3913 mul = MUL_NRM;
3914 }
3915 else if( nMatch( "TABL", chWord ) )
3916 {
3917 sd->sdCase = SD_TABLE;
3918 }
3919 else
3920 {
3921 sd->sdCase = SD_ILLEGAL;
3922 }
3923
3924 switch( sd->sdCase )
3925 {
3926 case SD_SINGLE_SIZE:
3927 /* single sized grain */
3928 mie_next_data(chFile,io2,chLine,&dl);
3929 mie_read_double(chFile,chLine,&sd->a[ipSize],true,dl);
3930 if( sd->a[ipSize] < SMALLEST_GRAIN || sd->a[ipSize] > LARGEST_GRAIN )
3931 {
3932 fprintf( ioQQQ, " Illegal value for grain size\n" );
3933 fprintf( ioQQQ, " Grain sizes should be between %.5f and %.0f micron\n",
3935 fprintf( ioQQQ, " Line #%ld: %s\n",dl,chLine);
3937 }
3938 break;
3939 case SD_NR_CARBON:
3940 /* single sized PAH with fixed number of carbon atoms */
3941 mie_next_data(chFile,io2,chLine,&dl);
3942 mie_read_long(chFile,chLine,&sd->nCarbon,true,dl);
3943 break;
3944 case SD_POWERLAW:
3945 /* simple power law distribution, first get lower limit */
3946 mie_next_data(chFile,io2,chLine,&dl);
3947 mie_read_double(chFile,chLine,&sd->a[ipBLo],true,dl);
3948 if( sd->a[ipBLo] < SMALLEST_GRAIN || sd->a[ipBLo] > LARGEST_GRAIN )
3949 {
3950 fprintf( ioQQQ, " Illegal value for grain size (lower limit)\n" );
3951 fprintf( ioQQQ, " Grain sizes should be between %.5f and %.0f micron\n",
3953 fprintf( ioQQQ, " Line #%ld: %s\n",dl,chLine);
3955 }
3956
3957 /* upper limit */
3958 mie_next_data(chFile,io2,chLine,&dl);
3959 mie_read_double(chFile,chLine,&sd->a[ipBHi],true,dl);
3960 if( sd->a[ipBHi] < SMALLEST_GRAIN || sd->a[ipBHi] > LARGEST_GRAIN ||
3961 sd->a[ipBHi] <= sd->a[ipBLo] )
3962 {
3963 fprintf( ioQQQ, " Illegal value for grain size (upper limit)\n" );
3964 fprintf( ioQQQ, " Grain sizes should be between %.5f and %.0f micron\n",
3966 fprintf( ioQQQ, " and upper limit should be greater than lower limit\n" );
3967 fprintf( ioQQQ, " Line #%ld: %s\n",dl,chLine);
3969 }
3970
3971 /* slope */
3972 mie_next_data(chFile,io2,chLine,&dl);
3973 if( sscanf( chLine, "%lf", &sd->a[ipExp] ) != 1 )
3974 {
3975 fprintf( ioQQQ, " Syntax error in %s\n",chFile);
3976 fprintf( ioQQQ, " Line #%ld: %s\n",dl,chLine);
3978 }
3979
3980 sd->a[ipBeta] = 0.;
3981 sd->a[ipSLo] = 0.;
3982 sd->a[ipSHi] = 0.;
3983
3984 sd->lim[ipBLo] = sd->a[ipBLo];
3985 sd->lim[ipBHi] = sd->a[ipBHi];
3986 break;
3987 case SD_EXP_CUTOFF1:
3988 case SD_EXP_CUTOFF2:
3989 case SD_EXP_CUTOFF3:
3990 /* powerlaw with first/second/third order exponential cutoff, inspired by
3991 * Greenberg (1978), Cosmic Dust, ed. J.A.M. McDonnell, Wiley, p. 187 */
3992 /* "lower limit", below this the exponential cutoff sets in */
3993 mie_next_data(chFile,io2,chLine,&dl);
3994 mie_read_double(chFile,chLine,&sd->a[ipBLo],false,dl);
3995
3996 /* "upper" limit, above this the exponential cutoff sets in */
3997 mie_next_data(chFile,io2,chLine,&dl);
3998 mie_read_double(chFile,chLine,&sd->a[ipBHi],false,dl);
3999
4000 /* exponent for power law */
4001 mie_next_data(chFile,io2,chLine,&dl);
4002 if( sscanf( chLine, "%lf", &sd->a[ipExp] ) != 1 )
4003 {
4004 fprintf( ioQQQ, " Syntax error in %s\n",chFile);
4005 fprintf( ioQQQ, " Line #%ld: %s\n",dl,chLine);
4007 }
4008
4009 /* beta parameter, for extra curvature in the powerlaw region */
4010 mie_next_data(chFile,io2,chLine,&dl);
4011 if( sscanf( chLine, "%lf", &sd->a[ipBeta] ) != 1 )
4012 {
4013 fprintf( ioQQQ, " Syntax error in %s\n",chFile);
4014 fprintf( ioQQQ, " Line #%ld: %s\n",dl,chLine);
4016 }
4017
4018 /* scale size for lower exponential cutoff, zero indicates normal cutoff */
4019 mie_next_data(chFile,io2,chLine,&dl);
4020 mie_read_double(chFile,chLine,&sd->a[ipSLo],false,dl);
4021
4022 /* scale size for upper exponential cutoff, zero indicates normal cutoff */
4023 mie_next_data(chFile,io2,chLine,&dl);
4024 mie_read_double(chFile,chLine,&sd->a[ipSHi],false,dl);
4025
4026 ref_neg = sd->a[ipBLo];
4027 step_neg = -mul*sd->a[ipSLo];
4028 ref_pos = sd->a[ipBHi];
4029 step_pos = mul*sd->a[ipSHi];
4030 lgTryOverride = true;
4031 break;
4032 case SD_LOG_NORMAL:
4033 /* log-normal distribution, first get center of gaussian */
4034 mie_next_data(chFile,io2,chLine,&dl);
4035 mie_read_double(chFile,chLine,&sd->a[ipGCen],true,dl);
4036
4037 /* 1-sigma width */
4038 mie_next_data(chFile,io2,chLine,&dl);
4039 mie_read_double(chFile,chLine,&sd->a[ipGSig],true,dl);
4040
4041 /* ref_pos, ref_neg is the grain radius at which a^4*dN/da peaks */
4042 ref_neg = ref_pos = sd->a[ipGCen]*exp(3.*pow2(sd->a[ipGSig]));
4043 step_neg = sd->a[ipGCen]*(exp(-mul*sd->a[ipGSig]) - 1.);
4044 step_pos = sd->a[ipGCen]*(exp(mul*sd->a[ipGSig]) - 1.);
4045 lgTryOverride = true;
4046 break;
4047 case SD_LIN_NORMAL:
4048 /* normal gaussian distribution, first get center of gaussian */
4049 mie_next_data(chFile,io2,chLine,&dl);
4050 mie_read_double(chFile,chLine,&sd->a[ipGCen],true,dl);
4051
4052 /* 1-sigma width */
4053 mie_next_data(chFile,io2,chLine,&dl);
4054 mie_read_double(chFile,chLine,&sd->a[ipGSig],true,dl);
4055
4056 /* ref_pos, ref_neg is the grain radius at which a^4*dN/da peaks */
4057 ref_neg = ref_pos = (sd->a[ipGCen] + sqrt(pow2(sd->a[ipGCen]) + 12.*pow2(sd->a[ipGSig])))/2.;
4058 step_neg = -mul*sd->a[ipGSig];
4059 step_pos = mul*sd->a[ipGSig];
4060 lgTryOverride = true;
4061 break;
4062 case SD_TABLE:
4063 /* user-supplied table of a^4*dN/da vs. a, first get lower limit on a */
4064 mie_next_data(chFile,io2,chLine,&dl);
4065 mie_read_double(chFile,chLine,&sd->a[ipBLo],true,dl);
4066 if( sd->a[ipBLo] < SMALLEST_GRAIN || sd->a[ipBLo] > LARGEST_GRAIN )
4067 {
4068 fprintf( ioQQQ, " Illegal value for grain size (lower limit)\n" );
4069 fprintf( ioQQQ, " Grain sizes should be between %.5f and %.0f micron\n",
4071 fprintf( ioQQQ, " Line #%ld: %s\n",dl,chLine);
4073 }
4074
4075 /* upper limit */
4076 mie_next_data(chFile,io2,chLine,&dl);
4077 mie_read_double(chFile,chLine,&sd->a[ipBHi],true,dl);
4078 if( sd->a[ipBHi] < SMALLEST_GRAIN || sd->a[ipBHi] > LARGEST_GRAIN ||
4079 sd->a[ipBHi] <= sd->a[ipBLo] )
4080 {
4081 fprintf( ioQQQ, " Illegal value for grain size (upper limit)\n" );
4082 fprintf( ioQQQ, " Grain sizes should be between %.5f and %.0f micron\n",
4084 fprintf( ioQQQ, " and upper limit should be greater than lower limit\n" );
4085 fprintf( ioQQQ, " Line #%ld: %s\n",dl,chLine);
4087 }
4088
4089 /* number of user supplied points */
4090 mie_next_data(chFile,io2,chLine,&dl);
4091 mie_read_long(chFile,chLine,&sd->npts,true,dl);
4092 if( sd->npts < 2 )
4093 {
4094 fprintf( ioQQQ, " Illegal value for no. of points in table\n" );
4095 fprintf( ioQQQ, " Line #%ld: %s\n",dl,chLine);
4097 }
4098
4099 /* allocate space for the table */
4100 sd->ln_a.resize(sd->npts);
4101 sd->ln_a4dNda.resize(sd->npts);
4102
4103 /* and read the table */
4104 for( i=0; i < sd->npts; ++i )
4105 {
4106 double help1, help2;
4107
4108 mie_next_data(chFile,io2,chLine,&dl);
4109 /* read data pair: a (micron), a^4*dN/da (arbitrary units) */
4110 if( sscanf( chLine, "%le %le", &help1, &help2 ) != 2 )
4111 {
4112 fprintf( ioQQQ, " Syntax error in %s\n",chFile);
4113 fprintf( ioQQQ, " Line #%ld: %s\n",dl,chLine);
4115 }
4116
4117 if( help1 <= 0. || help2 <= 0. )
4118 {
4119 fprintf( ioQQQ, " Reading table failed on line #%ld of %s\n",dl,chFile );
4120 fprintf( ioQQQ, " Illegal data value %.6e or %.6e\n", help1, help2 );
4122 }
4123
4124 sd->ln_a[i] = log(help1);
4125 sd->ln_a4dNda[i] = log(help2);
4126
4127 if( i > 0 && sd->ln_a[i] <= sd->ln_a[i-1] )
4128 {
4129 fprintf( ioQQQ, " Reading table failed on line #%ld of %s\n",dl,chFile );
4130 fprintf( ioQQQ, " Grain radii should be monotonically increasing\n" );
4132 }
4133 }
4134
4135 sd->lim[ipBLo] = sd->a[ipBLo];
4136 sd->lim[ipBHi] = sd->a[ipBHi];
4137 break;
4138 case SD_ILLEGAL:
4139 default:
4140 fprintf( ioQQQ, " unimplemented case for grain size distribution in file %s\n" , chFile );
4141 fprintf( ioQQQ, " Line #%ld: value %s\n",dl,chWord);
4143 }
4144
4145 /* >>chng 01 feb 12, use a^4*dN/da instead of dN/da to determine limits,
4146 * this assures that upper limit gives negligible mass fraction, PvH */
4147 /* in all cases where search_limit is used to determine lim[ipBLo] and lim[ipBHi],
4148 * the user can override these values in the last two lines of the size distribution
4149 * file. these inputs are mandatory, and should be given in the sequence lower
4150 * limit, upper limit. a value <= 0 indicates that search_limit should be used. */
4151 if( lgTryOverride )
4152 {
4153 double help;
4154
4155 mie_next_data(chFile,io2,chLine,&dl);
4156 mie_read_double(chFile,chLine,&help,false,dl);
4157 sd->lim[ipBLo] = ( help <= 0. ) ? search_limit(ref_neg,step_neg,FRAC_CUTOFF,*sd) : help;
4158
4159 mie_next_data(chFile,io2,chLine,&dl);
4160 mie_read_double(chFile,chLine,&help,false,dl);
4161 sd->lim[ipBHi] = ( help <= 0. ) ? search_limit(ref_pos,step_pos,FRAC_CUTOFF,*sd) : help;
4162
4163 if( sd->lim[ipBLo] < SMALLEST_GRAIN || sd->lim[ipBHi] > LARGEST_GRAIN ||
4164 sd->lim[ipBHi] <= sd->lim[ipBLo] )
4165 {
4166 fprintf( ioQQQ, " Illegal size limits: lower %.5f and/or upper %.5f\n",
4167 sd->lim[ipBLo], sd->lim[ipBHi] );
4168 fprintf( ioQQQ, " Grain sizes should be between %.5f and %.0f micron\n",
4170 fprintf( ioQQQ, " and upper limit should be greater than lower limit\n" );
4171 fprintf( ioQQQ, " Please alter the size distribution file\n" );
4173 }
4174 }
4175
4176 fclose(io2);
4177 return;
4178}
4179
4180
4181STATIC void mie_read_long(/*@in@*/ const char *chFile,
4182 /*@in@*/ const char chLine[], /* chLine[FILENAME_PATH_LENGTH_2] */
4183 /*@out@*/ long int *data,
4184 bool lgZeroIllegal,
4185 long int dl)
4186{
4187 DEBUG_ENTRY( "mie_read_long()" );
4188 if( sscanf( chLine, "%ld", data ) != 1 )
4189 {
4190 fprintf( ioQQQ, " Syntax error in %s\n",chFile);
4191 fprintf( ioQQQ, " Line #%ld: %s\n",dl,chLine);
4193 }
4194 if( *data < 0 || (*data == 0 && lgZeroIllegal) )
4195 {
4196 fprintf( ioQQQ, " Illegal data value in %s\n",chFile);
4197 fprintf( ioQQQ, " Line #%ld: %ld\n",dl,*data);
4199 }
4200 return;
4201}
4202
4203
4204STATIC void mie_read_realnum(/*@in@*/ const char *chFile,
4205 /*@in@*/ const char chLine[], /* chLine[FILENAME_PATH_LENGTH_2] */
4206 /*@out@*/ realnum *data,
4207 bool lgZeroIllegal,
4208 long int dl)
4209{
4210 DEBUG_ENTRY( "mie_read_realnum()" );
4211 double help;
4212 if( sscanf( chLine, "%lf", &help ) != 1 )
4213 {
4214 fprintf( ioQQQ, " Syntax error in %s\n",chFile);
4215 fprintf( ioQQQ, " Line #%ld: %s\n",dl,chLine);
4217 }
4218 *data = (realnum)help;
4219 if( *data < 0. || (*data == 0. && lgZeroIllegal) )
4220 {
4221 fprintf( ioQQQ, " Illegal data value in %s\n",chFile);
4222 fprintf( ioQQQ, " Line #%ld: %14.6e\n",dl,*data);
4224 }
4225 return;
4226}
4227
4228
4229STATIC void mie_read_double(/*@in@*/ const char *chFile,
4230 /*@in@*/ const char chLine[], /* chLine[FILENAME_PATH_LENGTH_2] */
4231 /*@out@*/ double *data,
4232 bool lgZeroIllegal,
4233 long int dl)
4234{
4235 DEBUG_ENTRY( "mie_read_double()" );
4236 if( sscanf( chLine, "%lf", data ) != 1 )
4237 {
4238 fprintf( ioQQQ, " Syntax error in %s\n",chFile);
4239 fprintf( ioQQQ, " Line #%ld: %s\n",dl,chLine);
4241 }
4242 if( *data < 0. || (*data == 0. && lgZeroIllegal) )
4243 {
4244 fprintf( ioQQQ, " Illegal data value in %s\n",chFile);
4245 fprintf( ioQQQ, " Line #%ld: %14.6e\n",dl,*data);
4247 }
4248 return;
4249}
4250
4251
4252STATIC void mie_read_form(/*@in@*/ const char chWord[], /* chWord[FILENAME_PATH_LENGTH_2] */
4253 /*@out@*/ double elmAbun[], /* elmAbun[LIMELM] */
4254 /*@out@*/ double *no_atoms,
4255 /*@out@*/ double *mol_weight)
4256{
4257 long int nelem,
4258 len;
4259 char chElmName[3];
4260 double frac;
4261
4262 DEBUG_ENTRY( "mie_read_form()" );
4263
4264 *no_atoms = 0.;
4265 *mol_weight = 0.;
4266 string strWord( chWord );
4267 for( nelem=0; nelem < LIMELM; nelem++ )
4268 {
4269 frac = 0.;
4270 strcpy(chElmName,elementnames.chElementSym[nelem]);
4271 if( chElmName[1] == ' ' )
4272 chElmName[1] = '\0';
4273 string::size_type ptr = strWord.find( chElmName );
4274 if( ptr != string::npos )
4275 {
4276 len = (long)strlen(chElmName);
4277 /* prevent spurious match, e.g. F matches Fe */
4278 if( !islower((unsigned char)chWord[ptr+len]) )
4279 {
4280 if( isdigit((unsigned char)chWord[ptr+len]) )
4281 {
4282 sscanf(chWord+ptr+len,"%lf",&frac);
4283 }
4284 else
4285 {
4286 frac = 1.;
4287 }
4288 }
4289 }
4290 elmAbun[nelem] = frac;
4291 /* >>chng 02 apr 22, don't count hydrogen in PAH's, PvH */
4292 if( nelem != ipHYDROGEN )
4293 *no_atoms += frac;
4294 *mol_weight += frac*dense.AtomicWeight[nelem];
4295 }
4296 /* prevent division by zero when no chemical formula was supplied */
4297 if( *no_atoms == 0. )
4298 *no_atoms = 1.;
4299 return;
4300}
4301
4302
4303STATIC void mie_write_form(/*@in@*/ const double elmAbun[], /* elmAbun[LIMELM] */
4304 /*@out@*/ char chWord[]) /* chWord[FILENAME_PATH_LENGTH_2] */
4305{
4306 long int len,
4307 nelem;
4308
4309 DEBUG_ENTRY( "mie_write_form()" );
4310
4311 len=0;
4312 for( nelem=0; nelem < LIMELM; nelem++ )
4313 {
4314 if( elmAbun[nelem] > 0. )
4315 {
4316 char chElmName[3];
4317 long index100 = nint(100.*elmAbun[nelem]);
4318
4319 strcpy(chElmName,elementnames.chElementSym[nelem]);
4320 if( chElmName[1] == ' ' )
4321 chElmName[1] = '\0';
4322
4323 if( index100 == 100 )
4324 sprintf( &chWord[len], "%s", chElmName );
4325 else if( index100%100 == 0 )
4326 sprintf( &chWord[len], "%s%ld", chElmName, index100/100 );
4327 else
4328 {
4329 double xIndex = (double)index100/100.;
4330 sprintf( &chWord[len], "%s%.2f", chElmName, xIndex );
4331 }
4332 len = (long)strlen( chWord );
4333 ASSERT( len < FILENAME_PATH_LENGTH_2-25 );
4334 }
4335 }
4336 return;
4337}
4338
4339
4340STATIC void mie_read_word(/*@in@*/ const char chLine[], /* chLine[FILENAME_PATH_LENGTH_2] */
4341 /*@out@*/ char chWord[], /* chWord[n] */
4342 long n,
4343 bool lgToUpper)
4344{
4345 long ip = 0,
4346 op = 0;
4347
4348 DEBUG_ENTRY( "mie_read_word()" );
4349
4350 /* skip leading spaces or double quotes */
4351 while( chLine[ip] == ' ' || chLine[ip] == '"' )
4352 ip++;
4353 /* now copy string until we hit next space or double quote */
4354 while( op < n-1 && chLine[ip] != ' ' && chLine[ip] != '"' )
4355 if( lgToUpper )
4356 chWord[op++] = toupper(chLine[ip++]);
4357 else
4358 chWord[op++] = chLine[ip++];
4359 /* the output string is always zero terminated */
4360 chWord[op] = '\0';
4361 return;
4362}
4363
4364/*=====================================================================*/
4365STATIC void mie_next_data(/*@in@*/ const char *chFile,
4366 /*@in@*/ FILE *io,
4367 /*@out@*/ char chLine[], /* chLine[FILENAME_PATH_LENGTH_2] */
4368 /*@in@*/ long int *dl)
4369{
4370 char *str;
4371
4372 DEBUG_ENTRY( "mie_next_data()" );
4373
4374 /* lines starting with a pound sign are considered comments and are skipped,
4375 * lines not starting with a pound sign are considered to contain useful data.
4376 * however, comments may still be appended to the line and will be erased. */
4377
4378 chLine[0] = '#';
4379 while( chLine[0] == '#' )
4380 {
4381 mie_next_line(chFile,io,chLine,dl);
4382 }
4383
4384 /* erase comment part of the line */
4385 str = strstr_s(chLine,"#");
4386 if( str != NULL )
4387 *str = '\0';
4388 return;
4389}
4390
4391
4392/*=====================================================================*/
4393STATIC void mie_next_line(/*@in@*/ const char *chFile,
4394 /*@in@*/ FILE *io,
4395 /*@out@*/ char chLine[], /* chLine[FILENAME_PATH_LENGTH_2] */
4396 /*@in@*/ long int *dl)
4397{
4398 DEBUG_ENTRY( "mie_next_line()" );
4399
4400 if( read_whole_line( chLine, FILENAME_PATH_LENGTH_2, io ) == NULL )
4401 {
4402 fprintf( ioQQQ, " Could not read from %s\n",chFile);
4403 if( feof(io) )
4404 fprintf( ioQQQ, " EOF reached\n");
4405 fprintf( ioQQQ, " This grain data file does not have the expected format.\n");
4407 }
4408 (*dl)++;
4409 return;
4410}
4411
4412/*=====================================================================*
4413 *
4414 * The routines gauss_init and gauss_legendre were derived from the
4415 * program cmieuvx.f.
4416 *
4417 * Written by: P.G. Martin (CITA), based on the code described in
4418 * >>refer grain physics Hansen, J. E., Travis, L. D. 1974, Space Sci. Rev., 16, 527
4419 *
4420 * The algorithm in gauss_legendre was modified by Peter van Hoof to
4421 * avoid FP overflow for large values of nn.
4422 *
4423 *=====================================================================*/
4424/* set up Gaussian quadrature for arbitrary interval */
4425void gauss_init(long int nn,
4426 double xbot,
4427 double xtop,
4428 const vector<double>& x, /* x[nn] */
4429 const vector<double>& a, /* a[nn] */
4430 vector<double>& rr, /* rr[nn] */
4431 vector<double>& ww) /* ww[nn] */
4432{
4433 long int i;
4434 double bma,
4435 bpa;
4436
4437 DEBUG_ENTRY( "gauss_init()" );
4438
4439 bpa = (xtop+xbot)/2.;
4440 bma = (xtop-xbot)/2.;
4441
4442 for( i=0; i < nn; i++ )
4443 {
4444 rr[i] = bpa + bma*x[nn-1-i];
4445 ww[i] = bma*a[i];
4446 }
4447 return;
4448}
4449
4450
4451/*=====================================================================*/
4452/* set up abscissas and weights for Gauss-Legendre intergration of arbitrary even order */
4453void gauss_legendre(long int nn,
4454 vector<double>& x, /* x[nn] */
4455 vector<double>& a) /* a[nn] */
4456{
4457 long int i,
4458 iter,
4459 j;
4460 double cc,
4461 csa,
4462 d,
4463 dp1,
4464 dpn = 0.,
4465 dq,
4466 fj,
4467 fn,
4468 pn,
4469 pn1 = 0.,
4470 q,
4471 xt = 0.;
4472
4473 const double SAFETY = 5.;
4474
4475
4476 DEBUG_ENTRY( "gauss_legendre()" );
4477
4478 if( nn%2 == 1 )
4479 {
4480 fprintf( ioQQQ, " Illegal number of abcissas\n" );
4482 }
4483
4484 vector<double> c(nn);
4485
4486 fn = (double)nn;
4487 csa = 0.;
4488 cc = 2.;
4489 for( j=1; j < nn; j++ )
4490 {
4491 fj = (double)j;
4492 /* >>chng 01 apr 10, prevent underflows in cc, pn, pn1, dpn and dp1 for large nn
4493 * renormalize c[j] -> 4*c[j], cc -> 4^(nn-1)*cc, hence cc = O(1), etc...
4494 * Old code: c[j] = pow2(fj)/(4.*(fj-0.5)*(fj+0.5)); */
4495 c[j] = pow2(fj)/((fj-0.5)*(fj+0.5));
4496 cc *= c[j];
4497 }
4498
4499 for( i=0; i < nn/2; i++ )
4500 {
4501 switch( i )
4502 {
4503 case 0:
4504 xt = 1. - 2.78/(4. + pow2(fn));
4505 break;
4506 case 1:
4507 xt = xt - 4.1*(1. + 0.06*(1. - 8./fn))*(1. - xt);
4508 break;
4509 case 2:
4510 xt = xt - 1.67*(1. + 0.22*(1. - 8./fn))*(x[0] - xt);
4511 break;
4512 default:
4513 xt = 3.*(x[i-1] - x[i-2]) + x[i-3];
4514 }
4515 d = 1.;
4516 for( iter=1; (iter < 20) && (fabs(d) > DBL_EPSILON); iter++ )
4517 {
4518 /* >>chng 01 apr 10, renormalize pn -> 2^(nn-1)*pn, dpn -> 2^(nn-1)*dpn
4519 * pn1 -> 2^(nn-2)*pn1, dp1 -> 2^(nn-2)*dp1
4520 * Old code: pn1 = 1.; */
4521 pn1 = 0.5;
4522 pn = xt;
4523 dp1 = 0.;
4524 dpn = 1.;
4525 for( j=1; j < nn; j++ )
4526 {
4527 /* >>chng 01 apr 10, renormalize pn -> 2^(nn-1)*pn, dpn -> 2^(nn-1)*dpn
4528 * Old code: q = xt*pn - c[j]*pn1; dq = xt*dpn - c[j]*dp1 + pn; */
4529 q = 2.*xt*pn - c[j]*pn1;
4530 dq = 2.*xt*dpn - c[j]*dp1 + 2.*pn;
4531 pn1 = pn;
4532 pn = q;
4533 dp1 = dpn;
4534 dpn = dq;
4535 }
4536 d = pn/dpn;
4537 xt -= d;
4538 }
4539 x[i] = xt;
4540 x[nn-1-i] = -xt;
4541 /* >>chng 01 apr 10, renormalize dpn -> 2^(nn-1)*dpn, pn1 -> 2^(nn-2)*pn1
4542 * Old code: a[i] = cc/(dpn*pn1); */
4543 a[i] = cc/(dpn*2.*pn1);
4544 a[nn-1-i] = a[i];
4545 csa += a[i];
4546 }
4547
4548 /* this routine has been tested for every even nn between 2 and 4096
4549 * it passed the test for each of those cases with SAFETY < 3.11 */
4550 if( fabs(1.-csa) > SAFETY*fn*DBL_EPSILON )
4551 {
4552 fprintf( ioQQQ, " gauss_legendre failed to converge: delta = %.4e\n", fabs(1.-csa) );
4554 }
4555 return;
4556}
4557
4558
4559/* find index ind such that min(xa[ind],xa[ind+1]) <= x <= max(xa[ind],xa[ind+1]).
4560 * xa is assumed to be strictly monotically increasing or decreasing.
4561 * if x is outside the range spanned by xa, lgOutOfBounds is raised and ind is set to -1
4562 * n is the number of elements in xa. */
4563void find_arr(double x,
4564 const vector<double>& xa,
4565 long int n,
4566 /*@out@*/ long int *ind,
4567 /*@out@*/ bool *lgOutOfBounds)
4568{
4569 long int i1,
4570 i2,
4571 i3,
4572 sgn,
4573 sgn2;
4574
4575 DEBUG_ENTRY( "find_arr()" );
4576 /* this routine works for strictly monotically increasing
4577 * and decreasing arrays, sgn indicates which case it is */
4578 if( n < 2 )
4579 {
4580 fprintf( ioQQQ, " Invalid array\n");
4582 }
4583
4584 i1 = 0;
4585 i3 = n-1;
4586 sgn = sign3(xa[i3]-xa[i1]);
4587 if( sgn == 0 )
4588 {
4589 fprintf( ioQQQ, " Ill-ordered array\n");
4591 }
4592
4593 *lgOutOfBounds = x < min(xa[0],xa[n-1]) || x > max(xa[0],xa[n-1]);
4594 if( *lgOutOfBounds )
4595 {
4596 *ind = -1;
4597 return;
4598 }
4599
4600 i2 = (n-1)/2;
4601 while( (i3-i1) > 1 )
4602 {
4603 sgn2 = sign3(x-xa[i2]);
4604 if( sgn2 != 0 )
4605 {
4606 if( sgn == sgn2 )
4607 {
4608 i1 = i2;
4609 }
4610 else
4611 {
4612 i3 = i2;
4613 }
4614 i2 = (i1+i3)/2;
4615 }
4616 else
4617 {
4618 *ind = i2;
4619 return;
4620 }
4621 }
4622 *ind = i1;
4623 return;
4624}
4625
4626
4627/*=====================================================================*
4628 *
4629 * The routines sinpar, anomal, bigk, ritodf, and dftori were derived
4630 * from the program cmieuvx.f.
4631 *
4632 * Written by: P.G. Martin (CITA), based on the code described in
4633 * >>refer grain physics Hansen, J. E., Travis, L. D. 1974, Space Sci. Rev., 16, 527
4634 *
4635 *=====================================================================*/
4636
4637/* Oct 1988 for UV - X-ray extinction, including anomalous diffraction check
4638 * this version reads in real and imaginary parts of the refractive
4639 * index, with imaginary part positive (nridf = 3) or nr-1 (nridf = 2) or
4640 * real and imaginary parts of the dielectric function (nridf = 1)
4641 * Dec 1988: added qback; approximation for small x;
4642 * qphase, better convergence checking
4643 *
4644 * in anomalous diffraction: qext and qabs calculated - qscatt by subtraction
4645 * in rayleigh-gans: qscatt and qabs calculated
4646 * in mie: qext and qscatt calculated
4647 * */
4648
4649/* sinpar.f
4650 * consistency checks updated july 1999
4651 * t1 updated mildly 19 oct 1992
4652 * utility for mieuvx.f and mieuvxsd.f */
4653static int NMXLIM = 16000;
4654
4655STATIC void sinpar(double nre,
4656 double nim,
4657 double x,
4658 /*@out@*/ double *qext,
4659 /*@out@*/ double *qphase,
4660 /*@out@*/ double *qscat,
4661 /*@out@*/ double *ctbrqs,
4662 /*@out@*/ double *qback,
4663 /*@out@*/ long int *iflag)
4664{
4665 long int n,
4666 nmx1,
4667 nmx2,
4668 nn,
4669 nsqbk;
4670 double ectb,
4671 eqext,
4672 eqpha,
4673 eqscat,
4674 error=0.,
4675 error1=0.,
4676 rx,
4677 t1,
4678 t2,
4679 t3,
4680 t4,
4681 t5,
4682 tx,
4683 x3,
4684 x5=0.,
4685 xcut,
4686 xrd;
4687 complex<double> cdum1,
4688 cdum2,
4689 ci,
4690 eqb,
4691 nc,
4692 nc2,
4693 nc212,
4694 qbck,
4695 rrf,
4696 rrfx,
4697 sman,
4698 sman1,
4699 smbn,
4700 smbn1,
4701 tc1,
4702 tc2,
4703 wn,
4704 wn1,
4705 wn2;
4706
4707 DEBUG_ENTRY( "sinpar()" );
4708
4709 vector< complex<double> > a(NMXLIM);
4710
4711 *iflag = 0;
4712 ci = complex<double>(0.,1.);
4713 nc = complex<double>(nre,-nim);
4714 nc2 = nc*nc;
4715 rrf = 1./nc;
4716 rx = 1./x;
4717 rrfx = rrf*rx;
4718
4719 /* t1 is the number of terms nmx2 that will be needed to obtain convergence
4720 * try to minimize this, because the a(n) downwards recursion has to
4721 * start at nmx1 larger than this
4722 *
4723 * major loop series is summed to nmx2, or less when converged
4724 * nmx1 is used for a(n) only, n up to nmx2.
4725 * must start evaluation sufficiently above nmx2 that a(nmx1)=(0.,0.)
4726 * is a good approximation
4727 *
4728 *
4729 *orig with slight modification for extreme UV and X-ray, n near 1., large x
4730 *orig t1=x*dmax1( 1.1d0,dsqrt(nr*nr+ni*ni) )*
4731 *orig 1(1.d0+0.02d0*dmax1(dexp(-x/100.d0)*x/10.d0,dlog10(x)))
4732 *
4733 * rules like those of wiscombe 1980 are slightly more efficient */
4734 xrd = exp(log(x)/3.);
4735 /* the final number in t1 was 1., 2. for large x, and 3. is needed sometimes
4736 * see also idnint use below */
4737 t1 = x + 4.*xrd + 3.;
4738 /* t1=t1+0.05d0*xrd
4739 * was 0., then 1., then 2., now 3. for intermediate x
4740 * 19 oct 1992 */
4741 if( !(x <= 8. || x >= 4200.) )
4742 t1 += 0.05*xrd + 3.;
4743 t1 *= 1.01;
4744
4745 /* the original rule of dave for starting the downwards recursion was
4746 * to start at 1.1*|mx| + 1, i.e. with the original version of t1
4747 *orig nmx1=1.10d0*t1
4748 *
4749 * try a simpler, less costly one, as in bohren and huffman, p 478
4750 * this is the form for use with wiscombe rules for t1
4751 * tests: it produces the same results as the more costly version
4752 * */
4753 t4 = x*sqrt(nre*nre+nim*nim);
4754 nmx1 = nint(max(t1,t4)) + 15;
4755
4756 if( nmx1 < NMXLIM )
4757 {
4758 nmx2 = nint(t1);
4759 /*orig if( nmx1 .gt. 150 ) go to 22
4760 *orig nmx1 = 150
4761 *orig nmx2 = 135
4762 *
4763 * try a more efficient scheme */
4764 if( nmx2 <= 4 )
4765 {
4766 nmx2 = 4;
4767 nmx1 = nint(max(4.,t4)) + 15;
4768 }
4769
4770 /* downwards recursion for logarithmic derivative */
4771 a[nmx1] = 0.;
4772
4773 /* note that with the method of lentz 1976 (appl opt 15, 668), it would be
4774 * possible to find a(nmx2) directly, and start the downwards recursion there
4775 * however, there is not much in it with above form for nmx1 which uses just */
4776 for( n=0; n < nmx1; n++ )
4777 {
4778 nn = nmx1 - n;
4779 a[nn-1] = (double)(nn+1)*rrfx - 1./((double)(nn+1)*rrfx+a[nn]);
4780 }
4781
4782 t1 = cos(x);
4783 t2 = sin(x);
4784 wn2 = complex<double>(t1,-t2);
4785 wn1 = complex<double>(t2,t1);
4786 wn = rx*wn1 - wn2;
4787 tc1 = a[0]*rrf + rx;
4788 tc2 = a[0]*nc + rx;
4789 sman = (tc1*wn.real() - wn1.real())/(tc1*wn - wn1);
4790 smbn = (tc2*wn.real() - wn1.real())/(tc2*wn - wn1);
4791
4792 /* small x; above calculations subject to rounding errors
4793 * see bohren and huffman p 131
4794 * wiscombe 1980 appl opt 19, 1505 gives alternative formulation */
4795 xcut = 3.e-04;
4796 if( x < xcut )
4797 {
4798 nc212 = (nc2-1.)/(nc2+2.);
4799 x3 = pow3(x);
4800 x5 = x3*pow2(x);
4801 /* note change sign convention for m = n - ik here */
4802 sman = ci*2.*x3*nc212*(1./3.+x*x*0.2*(nc2-2.)/(nc2+2.)) + 4.*x5*x*nc212*nc212/9.;
4803 smbn = ci*x5*(nc2-1.)/45.;
4804 }
4805
4806 sman1 = sman;
4807 smbn1 = smbn;
4808 t1 = 1.5;
4809 sman *= t1;
4810 smbn *= t1;
4811 /* case n=1; note previous multiplication of sman and smbn by t1=1.5 */
4812 *qext = 2.*(sman.real() + smbn.real());
4813 *qphase = 2.*(sman.imag() + smbn.imag());
4814 nsqbk = -1;
4815 qbck = -2.*(sman - smbn);
4816 *qscat = (norm(sman) + norm(smbn))/.75;
4817
4818 *ctbrqs = 0.0;
4819 n = 2;
4820
4821 /************************* Major loop begins here ************************/
4822 while( true )
4823 {
4824 t1 = 2.*(double)n - 1.;
4825 t3 = 2.*(double)n + 1.;
4826 wn2 = wn1;
4827 wn1 = wn;
4828 wn = t1*rx*wn1 - wn2;
4829 cdum1 = a[n-1];
4830 cdum2 = n*rx;
4831 tc1 = cdum1*rrf + cdum2;
4832 tc2 = cdum1*nc + cdum2;
4833 sman = (tc1*wn.real() - wn1.real())/(tc1*wn - wn1);
4834 smbn = (tc2*wn.real() - wn1.real())/(tc2*wn - wn1);
4835
4836 /* small x, n=2
4837 * see bohren and huffman p 131 */
4838 if( x < xcut && n == 2 )
4839 {
4840 /* note change sign convention for m = n - ik here */
4841 sman = ci*x5*(nc2-1.)/(15.*(2.*nc2+3.));
4842 smbn = 0.;
4843 }
4844
4845 eqext = t3*(sman.real() + smbn.real());
4846 *qext += eqext;
4847 eqpha = t3*(sman.imag() + smbn.imag());
4848 *qphase += eqpha;
4849 nsqbk = -nsqbk;
4850 eqb = t3*(sman - smbn)*(double)nsqbk;
4851 qbck += eqb;
4852 tx = norm(sman) + norm(smbn);
4853 eqscat = t3*tx;
4854 *qscat += eqscat;
4855 t2 = (double)(n - 1);
4856 t5 = (double)n;
4857 t4 = t1/(t5*t2);
4858 t2 = (t2*(t5 + 1.))/t5;
4859 ectb = t2*(sman1.real()*sman.real()+sman1.imag()*sman.imag() + smbn1.real()*smbn.real() +
4860 smbn1.imag()*smbn.imag()) +
4861 t4*(sman1.real()*smbn1.real()+sman1.imag()*smbn1.imag());
4862 *ctbrqs += ectb;
4863
4864 /* check convergence
4865 * could decrease for large x and small m-1 in UV - X-ray; probably negligible */
4866 if( tx < 1.e-14 )
4867 {
4868 /* looks good but check relative convergence */
4869 eqext = fabs(eqext/ *qext);
4870 eqpha = fabs(eqpha/ *qphase);
4871 eqscat = fabs(eqscat/ *qscat);
4872 ectb = ( n == 2 ) ? 0. : fabs(ectb/ *ctbrqs);
4873 eqb = complex<double>( fabs(eqb.real()/qbck.real()), fabs(eqb.imag()/qbck.imag()) );
4874 /* leave out eqb.re/im, which are sometimes least well converged */
4875 error = MAX4(eqext,eqpha,eqscat,ectb);
4876 /* put a milder constraint on eqb.re/im */
4877 error1 = max(eqb.real(),eqb.imag());
4878 if( error < 1.e-07 && error1 < 1.e-04 )
4879 break;
4880
4881 /* not sufficiently converged
4882 *
4883 * cut out after n=2 for small x, since approximation is being used */
4884 if( x < xcut )
4885 break;
4886 }
4887
4888 smbn1 = smbn;
4889 sman1 = sman;
4890 n++;
4891 if( n > nmx2 )
4892 {
4893 *iflag = 1;
4894 break;
4895 }
4896 }
4897 /* renormalize */
4898 t1 = 2.*pow2(rx);
4899 *qext *= t1;
4900 *qphase *= t1;
4901 *qback = norm(qbck)*pow2(rx);
4902 *qscat *= t1;
4903 *ctbrqs *= 2.*t1;
4904 }
4905 else
4906 {
4907 *iflag = 2;
4908 }
4909 return;
4910}
4911
4912
4913STATIC void anomal(double x,
4914 /*@out@*/ double *qext,
4915 /*@out@*/ double *qabs,
4916 /*@out@*/ double *qphase,
4917 /*@out@*/ double *xistar,
4918 double delta,
4919 double beta)
4920{
4921 /*
4922 *
4923 * in anomalous diffraction: qext and qabs calculated - qscatt by subtraction
4924 * in rayleigh-gans: qscatt and qabs calculated
4925 * in mie: qext and qscatt calculated
4926 *
4927 */
4928 double xi,
4929 xii;
4930 complex<double> cbigk,
4931 ci,
4932 cw;
4933
4934 DEBUG_ENTRY( "anomal()" );
4935 /* anomalous diffraction: x>>1 and |m-1|<<1, any xi,xii
4936 * original approach see Martin 1970. MN 149, 221 */
4937 xi = 2.*x*delta;
4938 xii = 2.*x*beta;
4939 /* xistar small is the basis for rayleigh-gans, any x, m-1 */
4940 *xistar = sqrt(pow2(xi)+pow2(xii));
4941 /* alternative approach see martin 1978 p 23 */
4942 ci = complex<double>(0.,1.);
4943 cw = -complex<double>(xi,xii)*ci;
4944 bigk(cw,&cbigk);
4945 *qext = 4.*cbigk.real();
4946 *qphase = 4.*cbigk.imag();
4947 cw = 2.*xii;
4948 bigk(cw,&cbigk);
4949 *qabs = 2.*cbigk.real();
4950 /* ?? put g in here - analytic version not known */
4951 return;
4952}
4953
4954
4955STATIC void bigk(complex<double> cw,
4956 /*@out@*/ complex<double> *cbigk)
4957{
4958 /*
4959 * see martin 1978 p 23
4960 */
4961
4962 DEBUG_ENTRY( "bigk()" );
4963 /* non-vax; use generic function */
4964 if( abs(cw) < 1.e-2 )
4965 {
4966 /* avoid severe loss of precision for small cw; expand exponential
4967 * coefficients are 1/n! - 1/(n+1)! = 1/(n+1)(n-1)!;n=2,3,4,5,6,7
4968 * accurate to (1+ order cw**6) */
4969 *cbigk = cw*((1./3.)-cw*((1./8.)-cw*((1./30.)-cw*((1./144.)-cw*((1./840.)-cw*(1./5760.))))));
4970 }
4971 else
4972 {
4973 *cbigk = 0.5 + (exp(-cw)*(1.+cw)-1.)/(cw*cw);
4974 }
4975 return;
4976}
4977
4978
4979/* utility for use with mieuvx/sd */
4980STATIC void ritodf(double nr,
4981 double ni,
4982 /*@out@*/ double *eps1,
4983 /*@out@*/ double *eps2)
4984{
4985
4986 DEBUG_ENTRY( "ritodf()" );
4987 /* refractive index to dielectric function */
4988 *eps1 = nr*nr - ni*ni;
4989 *eps2 = 2.*nr*ni;
4990 return;
4991}
4992
4993
4994/* utility for use with mieuvx/sd */
4995STATIC void dftori(/*@out@*/ double *nr,
4996 /*@out@*/ double *ni,
4997 double eps1,
4998 double eps2)
4999{
5000 double eps;
5001
5002 DEBUG_ENTRY( "dftori()" );
5003 /* dielectric function to refractive index */
5004 eps = sqrt(eps2*eps2+eps1*eps1);
5005 *nr = sqrt((eps+eps1)/2.);
5006 ASSERT( *nr > 0. );
5007 /* >>chng 03 jan 02, old expression for ni suffered
5008 * from cancellation error in the X-ray regime, PvH */
5009 /* *ni = sqrt((eps-eps1)/2.); */
5010 *ni = eps2/(2.*(*nr));
5011 return;
5012}
@ PHFIT95
Definition atmdat.h:276
static double a2[63]
static double x2[63]
static double x0[83]
static double x1[83]
static double a1[83]
t_called called
Definition called.cpp:5
FILE * ioQQQ
Definition cddefines.cpp:7
const int FILENAME_PATH_LENGTH_2
Definition cddefines.h:249
#define ASSERT(exp)
Definition cddefines.h:578
#define STATIC
Definition cddefines.h:97
long nMatch(const char *chKey, const char *chCard)
Definition service.cpp:451
const int ipCARBON
Definition cddefines.h:310
T pow2(T a)
Definition cddefines.h:931
int sign3(T x)
Definition cddefines.h:808
const int LIMELM
Definition cddefines.h:258
#define EXIT_FAILURE
Definition cddefines.h:140
char toupper(char c)
Definition cddefines.h:700
#define cdEXIT(FAIL)
Definition cddefines.h:434
const char * strstr_s(const char *haystack, const char *needle)
Definition cddefines.h:1429
float realnum
Definition cddefines.h:103
const char * strchr_s(const char *s, int c)
Definition cddefines.h:1439
long nint(double x)
Definition cddefines.h:719
#define POW4
Definition cddefines.h:943
long max(int a, long b)
Definition cddefines.h:775
bool fp_equal(sys_float x, sys_float y, int n=3)
Definition cddefines.h:812
char * read_whole_line(char *chLine, int nChar, FILE *ioIN)
Definition service.cpp:70
NORETURN void TotalInsanity(void)
Definition service.cpp:886
#define MAX4(a, b, c, d)
Definition cddefines.h:792
const int ipHYDROGEN
Definition cddefines.h:305
#define DEBUG_ENTRY(funcname)
Definition cddefines.h:684
long min(int a, long b)
Definition cddefines.h:723
T pow3(T a)
Definition cddefines.h:938
void ShowMe(void)
Definition service.cpp:181
double powi(double, long int)
Definition service.cpp:604
static t_version & Inst()
Definition cddefines.h:175
vector< double > wavlen[NAX]
vector< double > nr1[NAX]
mat_type matType
double wt[NAX]
long int charge
double atom_weight
long int nOpcData
double mol_weight
vector< complex< double > > n[NAX]
double subl_temp
double therm_eff
long int ndata[NAX]
long int cAxis
long int magic
double elmAbun[LIMELM]
void p_clear0()
vector< double > opcAnu
vector< double > opcData[NDAT]
double bandgap
void p_clear1()
long int nOpcCols
rfi_type rfiType
long int nAxes
const multi_geom< d, ALLOC > & clone() const
double area
long int nCarbon
vector< double > ln_a
double unity
long int nn
long int nmul
sd_type sdCase
vector< double > xx
void p_clear1()
double radius
double unity_bin
vector< double > rr
double cSize
long int nPart
double a[NSD]
vector< double > aa
long int npts
long int cPart
void clear()
vector< double > ww
double clim[2]
bool lgLogScale
vector< double > ln_a4dNda
double lim[2]
double vol
long int magic
double phfit(long int nz, long int ne, long int is, double e)
void set_version(phfit_version val)
Definition atmdat.h:318
t_continuum continuum
Definition continuum.cpp:5
FILE * open_data(const char *fname, const char *mode, access_scheme scheme)
Definition cpu.cpp:625
static t_cpu cpu
Definition cpu.h:355
@ AS_LOCAL_ONLY
Definition cpu.h:208
@ AS_DATA_LOCAL
Definition cpu.h:208
t_dense dense
Definition dense.cpp:24
t_elementnames elementnames
static const double TOLER
Definition grains.cpp:77
static const int LABELSUB2
void mie_read_opc(const char *chFile, const GrainPar &gp)
static int NMXLIM
static const long MAGIC_OPC
STATIC void init_eps(double, long, const vector< grain_data > &, vector< complex< double > > &)
static const bool pah3_hoc[30]
STATIC void anomal(double, double *, double *, double *, double *, double, double)
STATIC void dftori(double *, double *, double, double)
STATIC void bigk(complex< double >, complex< double > *)
STATIC void mie_read_long(const char *, const char[], long int *, bool, long int)
static const int NDAT
static const int ipAlpha
STATIC void mie_read_realnum(const char *, const char[], realnum *, bool, long int)
STATIC void mie_integrate(sd_data *, double, double, double *)
static const int ipBeta
void car3_fun(double, const sd_data *, const grain_data[], double *, double *, double *, int *)
static const double pah3_width[30]
static const double pah3n_strength[30]
void mie_write_opc(const char *rfi_file, const char *szd_file, long int nbin)
STATIC void Stognienko(complex< double >, const vector< double > &, const vector< complex< double > > &, long, complex< double > *, double *, double *)
static const long MAGIC_SZD
sd_type
@ SD_LIN_NORMAL
@ SD_LOG_NORMAL
@ SD_EXP_CUTOFF2
@ SD_EXP_CUTOFF3
@ SD_NR_CARBON
@ SD_ILLEGAL
@ SD_TABLE
@ SD_SINGLE_SIZE
@ SD_POWERLAW
@ SD_EXP_CUTOFF1
STATIC void mie_next_line(const char *, FILE *, char *, long int *)
static const int LABELSIZE
emt_type
@ FARAFONOV00
@ BRUGGEMAN35
@ STOGNIENKO95
static const double LARGEST_GRAIN
static const int NAX
STATIC void mie_write_form(const double[], char[])
STATIC complex< double > cnewton(void(*)(complex< double >, const vector< double > &, const vector< complex< double > > &, long, complex< double > *, double *, double *), const vector< double > &, const vector< complex< double > > &, long, complex< double >, double)
static const double pah1_width[7]
STATIC void pah3_fun(double, const sd_data *, const grain_data *, double *, double *, double *, int *)
void gauss_legendre(long int nn, vector< double > &x, vector< double > &a)
static const int ipSHi
static const double pah2n_strength[14]
static const int LABELSUB1
STATIC void tbl_fun(double, const sd_data *, const grain_data *, double *, double *, double *, int *)
STATIC void ritodf(double, double, double *, double *)
STATIC void mie_next_data(const char *, FILE *, char *, long int *)
STATIC void mie_read_word(const char[], char[], long, bool)
static const double pah2_wavl[14]
STATIC void mie_read_szd(const char *, sd_data *)
STATIC void pah2_fun(double, const sd_data *, const grain_data *, double *, double *, double *, int *)
static const int ipGSig
const int NPTS_DERIV
STATIC void mie_repair(const char *, long, int, int, const double[], double[], vector< int > &, bool, bool *)
STATIC double mie_find_slope(const double[], const double[], const vector< int > &, long, long, int, bool, bool *)
double Drude(double, double, double, double)
static const double pah3_wavl[30]
static const int ipExp
static const double SMALLEST_GRAIN
static const int ipBLo
static const long MAGIC_MIX
STATIC void pah1_fun(double, const sd_data *, const grain_data *, double *, double *, double *, int *)
STATIC void mie_cs_size_distr(double, sd_data *, const grain_data *, void(*)(double, const sd_data *, const grain_data *, double *, double *, double *, int *), double *, double *, double *, int *)
STATIC void mie_read_form(const char *, double[], double *, double *)
static const long MIX_TABLE_SIZE
rfi_type
@ OPC_PAH3N
@ OPC_PAH2C
@ OPC_PAH1
@ OPC_PAH2N
@ OPC_TABLE
@ RFI_TABLE
@ OPC_GREY
@ OPC_PAH3C
static const double pah2_width[14]
static const double pah1_strength[7]
static const int NSD
void car2_fun(double, const sd_data *, const grain_data[], double *, double *, double *, int *)
STATIC void mie_read_mix(const char *, grain_data *)
STATIC void mie_read_rfi(const char *, grain_data *)
void find_arr(double x, const vector< double > &xa, long int n, long int *ind, bool *lgOutOfBounds)
STATIC double search_limit(double, double, double, sd_data)
STATIC void sinpar(double, double, double, double *, double *, double *, double *, double *, long *)
STATIC void mie_calc_ial(const grain_data *, long, vector< double > &, const char *, bool *)
void car1_fun(double, const sd_data *, const grain_data[], double *, double *, double *, int *)
static const double pah3c_strength[30]
STATIC void ld01_fun(void(*)(double, const sd_data *, const grain_data[], double *, double *, double *, int *), double, double, double, const sd_data *, const grain_data[], double *, double *, double *, int *)
static const int ipSize
static const int ipGCen
static const int ipSLo
const int NPTS_COMB
STATIC void mie_auxiliary(sd_data *, const grain_data *, const char *)
static const int ipBHi
static const long MAGIC_RFI
STATIC void mie_read_double(const char *, const char[], double *, bool, long int)
void gauss_init(long int nn, double xbot, double xtop, const vector< double > &x, const vector< double > &a, vector< double > &rr, vector< double > &ww)
STATIC double size_distr(double, const sd_data *)
static const double pah2c_strength[14]
STATIC bool mie_auxiliary2(vector< int > &, multi_arr< double, 2 > &, multi_arr< double, 2 > &, multi_arr< double, 2 > &, long, long)
static const int WORDLEN
STATIC void Bruggeman(complex< double >, const vector< double > &, const vector< complex< double > > &, long, complex< double > *, double *, double *)
STATIC void mie_cs(double, const sd_data *, const grain_data *, double *, double *, double *, int *)
static const double pah1_wlBand[7]
static const long LOOP_MAX
static const double SAFETY
double no_atoms(size_t nd)
GrainVar gv
Definition grainvar.cpp:5
mat_type
Definition grainvar.h:99
@ MAT_TOP
Definition grainvar.h:110
ial_type
Definition grainvar.h:67
@ IAL_SIL
Definition grainvar.h:69
@ IAL_CAR
Definition grainvar.h:68
diatomics h2("h2", 4100., &hmi.H2_total, Yan_H2_CS)
UNUSED const double PI
Definition physconst.h:29
UNUSED const double PI4
Definition physconst.h:35
UNUSED const double RYD_INF
Definition physconst.h:115
UNUSED const double EVRYD
Definition physconst.h:189
UNUSED const double ELECTRIC_CONST
Definition physconst.h:150
UNUSED const double ATOMIC_MASS_UNIT
Definition physconst.h:88
UNUSED const double WAVNRYD
Definition physconst.h:173
t_radius radius
Definition radius.cpp:5
t_rfield rfield
Definition rfield.cpp:8
static double * g
Definition species2.cpp:28
bool lgGreyGrain
Definition grainvar.h:117
bool lgRequestQHeating
Definition grainvar.h:118
df_type nDustFunc
Definition grainvar.h:115
bool lgForbidQHeating
Definition grainvar.h:116
double dep
Definition grainvar.h:114
static const int M
static const unsigned int NMD5
Definition thirdparty.h:384