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