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