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

Generated for cloudy by doxygen 1.8.1.1