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