00001
00002
00003 #include <ctype.h>
00004 #include "cddefines.h"
00005 #include "elementnames.h"
00006 #include "physconst.h"
00007 #include "path.h"
00008 #include "dense.h"
00009 #include "called.h"
00010 #include "version.h"
00011 #include "grainvar.h"
00012 #include "rfield.h"
00013 #include "atmdat.h"
00014 #include "grains.h"
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038 static const long MAGIC_RFI = 1030103L;
00039 static const long MAGIC_SZD = 2010403L;
00040 static const long MAGIC_OPC = 3030307L;
00041 static const long MAGIC_MIX = 4030103L;
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051 static const double SMALLEST_GRAIN = 0.0001*(1.-10.*DBL_EPSILON);
00052 static const double LARGEST_GRAIN = 10.*(1.+10.*DBL_EPSILON);
00053
00054
00055 static const int NSD = 7;
00056
00057
00058
00059 static const int ipSize = 0;
00060 static const int ipBLo = 0;
00061 static const int ipBHi = 1;
00062 static const int ipExp = 2;
00063 static const int ipBeta = 3;
00064 static const int ipSLo = 4;
00065 static const int ipSHi = 5;
00066 static const int ipAlpha = 6;
00067 static const int ipGCen = 2;
00068 static const int ipGSig = 3;
00069
00070
00071 typedef enum {
00072 RFI_TABLE, OPC_TABLE, OPC_GREY, OPC_PAH1
00073 } rfi_type;
00074
00075
00076 typedef enum {
00077 FARAFONOV00, STOGNIENKO95, BRUGGEMAN35
00078 } emt_type;
00079
00080
00081 typedef enum {
00082 SD_ILLEGAL, SD_SINGLE_SIZE, SD_POWERLAW, SD_EXP_CUTOFF1, SD_EXP_CUTOFF2,
00083 SD_EXP_CUTOFF3, SD_LOG_NORMAL, SD_LIN_NORMAL, SD_TABLE
00084 } sd_type;
00085
00086 typedef struct {
00087 double a[NSD];
00088 double lim[2];
00089 double clim[2];
00090 double *xx;
00091 double *aa;
00092 double *rr;
00093 double *ww;
00094 double unity;
00095 double unity_bin;
00096 double cSize;
00097 double radius;
00098 double area;
00099 double vol;
00100 double *ln_a;
00101 double *ln_a4dNda;
00102 sd_type sdCase;
00103 long int magic;
00104 long int cPart;
00105 long int nPart;
00106 long int nmul;
00107 long int nn;
00108 long int npts;
00109 bool lgLogScale;
00110 } sd_data;
00111
00112
00113 static const int NAX = 3;
00114 static const int NDAT = 4;
00115
00116 typedef struct {
00117 double *wavlen[NAX];
00118 complex<double> *n[NAX];
00119 double *nr1[NAX];
00120 double *opcAnu;
00121 double *opcData[NDAT];
00122 double wt[NAX];
00123 double abun;
00124 double depl;
00125 double elmAbun[LIMELM];
00126 double mol_weight;
00127 double atom_weight;
00128 double rho;
00129 double norm;
00130 double work;
00131 double bandgap;
00132 double therm_eff;
00133 double subl_temp;
00134 long int magic;
00135 long int cAxis;
00136 long int nAxes;
00137 long int ndata[NAX];
00138 long int nOpcCols;
00139 long int nOpcData;
00140 mat_type matType;
00141 rfi_type rfiType;
00142 } grain_data;
00143
00144
00145
00146
00147
00148 static const int WORDLEN = 5;
00149
00150
00151 static const int LABELSUB1 = 3;
00152 static const int LABELSUB2 = 5;
00153 static const int LABELSIZE = LABELSUB1 + LABELSUB2 + 4;
00154
00155
00156 static const long MIX_TABLE_SIZE = 2000L;
00157
00158 static void mie_auxiliary(sd_data*,const char*);
00159 static void mie_integrate(sd_data*,double,double,double*,bool);
00160 static void mie_cs_size_distr(double,sd_data*,grain_data*,
00161 void(*)(double,sd_data*,grain_data*,double*,
00162 double*,double*,int*),
00163 double*,double*,double*,int*);
00164 static void mie_cs(double,sd_data*,grain_data*,double*,double*,
00165 double*,int*);
00166 static void pah1_fun(double,sd_data*,grain_data*,double*,double*,
00167 double*,int*);
00168 static void tbl_fun(double,sd_data*,grain_data*,double*,double*,
00169 double*,int*);
00170 static double size_distr(double,sd_data*);
00171 static double search_limit(double,double,double,sd_data);
00172 static void mie_calc_ial(grain_data*,long,double[],const char*,bool*);
00173 static void mie_repair(const char*,long,int,int,float[],double[],int[],
00174 bool,bool*);
00175 static double mie_find_slope(const float[],const double[],const int[],
00176 long,long,int,bool,bool*);
00177 static void mie_read_rfi(const char*,grain_data*);
00178 static void mie_read_mix(const char*,grain_data*);
00179 static void init_eps(double,long,grain_data[],complex<double>[]);
00180 static complex<double> cnewton(
00181 void(*)(complex<double>,double[],complex<double>[],long,complex<double>*,double*,double*),
00182 double[],complex<double>[],long,complex<double>,double);
00183 static void Stognienko(complex<double>,double[],complex<double>[],long,complex<double>*,double*,double*);
00184 static void Bruggeman(complex<double>,double[],complex<double>[],long,complex<double>*,double*,double*);
00185 static void mie_read_szd(const char*,sd_data*);
00186 static void mie_read_long(const char*,const char[],long int*,bool,long int);
00187 static void mie_read_float(const char*,const char[],float*,bool,long int);
00188 static void mie_read_double(const char*,const char[],double*,bool,long int);
00189 static void mie_read_form(const char*,double[],double*,double*);
00190 static void mie_write_form(const double[],char[]);
00191 static void mie_read_word(const char[],char[],long,int);
00192 static void mie_next_data(const char*,FILE*,char*,long int*);
00193 static void mie_next_line(const char*,FILE*,char*,long int*);
00194
00195
00196
00197
00198 static void sinpar(double,double,double,double*,double*,double*,
00199 double*,double*,long*);
00200 static void anomal(double,double*,double*,double*,double*,double,double);
00201 static void bigk(complex<double>,complex<double>*);
00202 static void ritodf(double,double,double*,double*);
00203 static void dftori(double*,double*,double,double);
00204
00205
00206
00207
00208
00209 void mie_write_opc( const char *rfi_file,
00210 const char *szd_file,
00211 long int nbin)
00212 {
00213 int Error = 0,
00214
00215 *ErrorIndex;
00216 bool lgErr,
00217 lgErrorOccurred,
00218 lgWarning;
00219 long int i,
00220 j,
00221 nelem,
00222 p;
00223 double **acs_abs,
00224 **acs_sct,
00225 **a1g,
00226 *inv_att_len,
00227 cosb,
00228 cs_abs,
00229 cs_sct,
00230 volfrac,
00231 volnorm,
00232 wavlen;
00233 char chGrainLabel[LABELSIZE+1],
00234 ext[3],
00235 chFile[FILENAME_PATH_LENGTH_2],
00236 chFile2[FILENAME_PATH_LENGTH_2],
00237 hlp1[LABELSUB1+2],
00238 hlp2[LABELSUB2+2],
00239 *str,
00240 chString[FILENAME_PATH_LENGTH_2];
00241 sd_data sd;
00242 grain_data gd,
00243 gd2;
00244 time_t timer;
00245 FILE *fdes;
00246
00247
00248
00249 const long NPTS_TABLE = 100L;
00250
00251 DEBUG_ENTRY( "mie_write_opc()" );
00252
00253
00254 ErrorIndex = (int*)MALLOC(sizeof(int)*(unsigned)rfield.nupper);
00255
00256 gd.nAxes = 0;
00257 for( j=0; j < NAX; j++ )
00258 {
00259 gd.wavlen[j] = NULL;
00260 gd.n[j] = NULL;
00261 gd.nr1[j] = NULL;
00262 }
00263 gd.nOpcCols = 0;
00264 gd.opcAnu = NULL;
00265 for( j=0; j < NDAT; j++ )
00266 {
00267 gd.opcData[j] = NULL;
00268 }
00269 gd2.nAxes = 0;
00270 for( j=0; j < NAX; j++ )
00271 {
00272 gd2.wavlen[j] = NULL;
00273 gd2.n[j] = NULL;
00274 gd2.nr1[j] = NULL;
00275 }
00276 gd2.nOpcCols = 0;
00277 gd2.opcAnu = NULL;
00278 for( j=0; j < NDAT; j++ )
00279 {
00280 gd2.opcData[j] = NULL;
00281 }
00282 sd.ln_a = NULL;
00283 sd.ln_a4dNda = NULL;
00284
00285 mie_read_szd( szd_file , &sd );
00286
00287 sd.nPart = ( sd.sdCase == SD_SINGLE_SIZE ) ? 1 : nbin;
00288 if( sd.nPart <= 0 || sd.nPart >= 100 )
00289 {
00290 fprintf( ioQQQ, " Illegal number of size distribution bins: %ld\n",sd.nPart );
00291 fprintf( ioQQQ, " The number should be between 1 and 99.\n" );
00292 puts( "[Stop in mie_write_opc]" );
00293 cdEXIT(EXIT_FAILURE);
00294 }
00295 sd.lgLogScale = true;
00296
00297 string rfi_string ( rfi_file );
00298 if( rfi_string.find( ".rfi" ) != string::npos )
00299 {
00300 mie_read_rfi( rfi_file , &gd );
00301 }
00302 else if( rfi_string.find( ".mix" ) != string::npos )
00303 {
00304 mie_read_mix( rfi_file , &gd );
00305 }
00306 else
00307 {
00308 fprintf( ioQQQ, " Refractive index file name %s has wrong extention\n", rfi_file );
00309 fprintf( ioQQQ, " It should have extention .rfi or .mix.\n" );
00310 puts( "[Stop in mie_write_opc]" );
00311 cdEXIT(EXIT_FAILURE);
00312 }
00313
00314 if( gd.rfiType == OPC_TABLE && sd.nPart > 1 )
00315 {
00316 fprintf( ioQQQ, " Illegal number of size distribution bins: %ld\n",sd.nPart );
00317 fprintf( ioQQQ, " The number should always be 1 for OPC_TABLE files.\n" );
00318 puts( "[Stop in mie_write_opc]" );
00319 cdEXIT(EXIT_FAILURE);
00320 }
00321 if( gd.rho <= 0. )
00322 {
00323 fprintf( ioQQQ, " Illegal value for the specific density: %.4e\n",gd.rho );
00324 puts( "[Stop in mie_write_opc]" );
00325 cdEXIT(EXIT_FAILURE);
00326 }
00327 if( gd.mol_weight <= 0. )
00328 {
00329 fprintf( ioQQQ, " Illegal value for the molecular weight: %.4e\n",gd.mol_weight );
00330 puts( "[Stop in mie_write_opc]" );
00331 cdEXIT(EXIT_FAILURE);
00332 }
00333
00334 lgWarning = false;
00335
00336
00337 strcpy(chFile,rfi_file);
00338 str = strstr(chFile,".");
00339
00340 if( str != NULL )
00341 *str = '\0';
00342
00343 strcpy(chFile2,szd_file);
00344 str = strstr(chFile2,".");
00345
00346 if( str != NULL )
00347 *str = '\0';
00348
00349 if( sd.sdCase != SD_SINGLE_SIZE )
00350 {
00351 sprintf(ext,"%02ld",nbin);
00352 strcat(strcat(strcat(strcat(strcat(chFile,"_"),chFile2),"_"),ext),".opc");
00353 }
00354 else
00355 {
00356 strcat(strcat(strcat(chFile,"_"),chFile2),".opc");
00357 }
00358
00359 fprintf( ioQQQ, "\n Starting mie_write_opc, output will be written to %s\n\n",chFile );
00360
00361 mie_auxiliary(&sd,"init");
00362
00363
00364 gd.norm = sd.vol*gd.rho/(ATOMIC_MASS_UNIT*gd.mol_weight*gd.abun*gd.depl);
00365 volnorm = sd.vol;
00366 volfrac = 1.;
00367
00368 acs_abs = (double **)MALLOC(sizeof(double *)*(unsigned)sd.nPart);
00369 acs_sct = (double **)MALLOC(sizeof(double *)*(unsigned)sd.nPart);
00370 a1g = (double **)MALLOC(sizeof(double *)*(unsigned)sd.nPart);
00371 inv_att_len = (double *)MALLOC(sizeof(double)*(unsigned)rfield.nupper);
00372
00373 if( acs_abs == NULL || acs_sct == NULL || a1g == NULL || inv_att_len == NULL )
00374 {
00375 fprintf( ioQQQ, " Could not MALLOC arrays\n" );
00376 puts( "[Stop in mie_write_opc]" );
00377 cdEXIT(EXIT_FAILURE);
00378 }
00379 for( p=0; p < sd.nPart; p++ )
00380 {
00381 acs_abs[p] = (double *)MALLOC(sizeof(double)*(unsigned)rfield.nupper);
00382 acs_sct[p] = (double *)MALLOC(sizeof(double)*(unsigned)rfield.nupper);
00383 a1g[p] = (double *)MALLOC(sizeof(double)*(unsigned)rfield.nupper);
00384 if( acs_abs[p] == NULL || acs_sct[p] == NULL || a1g[p] == NULL )
00385 {
00386 fprintf( ioQQQ, " Could not MALLOC arrays\n" );
00387 puts( "[Stop in mie_write_opc]" );
00388 cdEXIT(EXIT_FAILURE);
00389 }
00390 }
00391
00392 if( (fdes = fopen( chFile, "w" )) == NULL )
00393 {
00394 fprintf( ioQQQ, " Could not open %s for writing\n", chFile );
00395 puts( "[Stop in mie_write_opc]" );
00396 cdEXIT(EXIT_FAILURE);
00397 }
00398 lgErr = false;
00399
00400 (void)time(&timer);
00401 lgErr = lgErr || ( fprintf(fdes,"# this file was created by Cloudy %s (%s) on %s",
00402 version.chVersion,version.chDate,ctime(&timer)) < 0 );
00403 lgErr = lgErr || ( fprintf(fdes,"# ===========================================\n#\n") < 0 );
00404 lgErr = lgErr || ( fprintf(fdes,"%12ld # magic number opacity file\n",MAGIC_OPC) < 0 );
00405 lgErr = lgErr || ( fprintf(fdes,"%12ld # magic number rfi/mix file\n",gd.magic) < 0 );
00406 lgErr = lgErr || ( fprintf(fdes,"%12ld # magic number szd file\n",sd.magic) < 0 );
00407
00408
00409
00410 strncpy(hlp1,chFile,(size_t)(LABELSUB1+1));
00411 hlp1[LABELSUB1+1] = '\0';
00412 str = strstr(hlp1,"-");
00413
00414 if( str != NULL )
00415 *str = '\0';
00416
00417 strncpy(hlp2,chFile2,(size_t)(LABELSUB2+1));
00418 hlp2[LABELSUB2+1] = '\0';
00419 str = strstr(hlp2,"-");
00420
00421 if( str != NULL )
00422 *str = '\0';
00423
00424 strcpy(chGrainLabel," ");
00425 if( sd.nPart > 1 )
00426 {
00427 hlp1[LABELSUB1] = '\0';
00428 hlp2[LABELSUB2] = '\0';
00429 strcat(strcat(strcat(strcat(chGrainLabel,hlp1),"-"),hlp2),"xx");
00430 lgErr = lgErr || ( fprintf(fdes,"%-12.12s # grain type label, xx will be replaced by bin no.\n",
00431 chGrainLabel) < 0 );
00432 }
00433 else
00434 {
00435 strcat(strcat(strcat(chGrainLabel,hlp1),"-"),hlp2);
00436 lgErr = lgErr || ( fprintf(fdes,"%-12.12s # grain type label\n", chGrainLabel) < 0 );
00437 }
00438
00439 lgErr = lgErr || ( fprintf(fdes,"%.6e # specific weight (g/cm^3)\n",gd.rho) < 0 );
00440 lgErr = lgErr || ( fprintf(fdes,"%.6e # molecular weight of grain molecule (amu)\n",gd.mol_weight) < 0 );
00441 lgErr = lgErr || ( fprintf(fdes,"%.6e # average molecular weight per atom (amu)\n", gd.atom_weight) < 0 );
00442 lgErr = lgErr || ( fprintf(fdes,"%.6e # abundance of grain molecule at max depletion\n",gd.abun) < 0 );
00443 lgErr = lgErr || ( fprintf(fdes,"%.6e # depletion efficiency\n",gd.depl) < 0 );
00444 lgErr = lgErr || ( fprintf(fdes,"%.6e # average grain radius <a^3>/<a^2>, full size distr (cm)\n",
00445 3.*sd.vol/sd.area) < 0 );
00446 lgErr = lgErr || ( fprintf(fdes,"%.6e # average grain surface area <4pi*a^2>, full size distr (cm^2)\n",
00447 sd.area) < 0 );
00448 lgErr = lgErr || ( fprintf(fdes,"%.6e # average grain volume <4/3pi*a^3>, full size distr (cm^3)\n",
00449 sd.vol) < 0 );
00450 lgErr = lgErr || ( fprintf(fdes,"%.6e # total grain radius Int(a) per H, full size distr (cm/H)\n",
00451 sd.radius/gd.norm) < 0 );
00452 lgErr = lgErr || ( fprintf(fdes,"%.6e # total grain area Int(4pi*a^2) per H, full size distr (cm^2/H)\n",
00453 sd.area/gd.norm) < 0 );
00454 lgErr = lgErr || ( fprintf(fdes,"%.6e # total grain vol Int(4/3pi*a^3) per H, full size distr (cm^3/H)\n",
00455 sd.vol/gd.norm) < 0 );
00456 lgErr = lgErr || ( fprintf(fdes,"%.6e # work function (Ryd)\n",gd.work) < 0 );
00457 lgErr = lgErr || ( fprintf(fdes,"%.6e # gap between valence and conduction band (Ryd)\n",gd.bandgap) < 0 );
00458 lgErr = lgErr || ( fprintf(fdes,"%.6e # efficiency of thermionic emission\n",gd.therm_eff) < 0 );
00459 lgErr = lgErr || ( fprintf(fdes,"%.6e # sublimation temperature (K)\n",gd.subl_temp) < 0 );
00460 lgErr = lgErr || ( fprintf(fdes,"%12d # material type, 1=carbonaceous, 2=silicate\n",gd.matType) < 0 );
00461 lgErr = lgErr || ( fprintf(fdes,"#\n# abundances of constituent elements rel. to hydrogen\n#\n") < 0 );
00462
00463 for( nelem=0; nelem < LIMELM; nelem++ )
00464 {
00465 lgErr = lgErr || ( fprintf(fdes,"%.6e # %s\n",gd.elmAbun[nelem],
00466 elementnames.chElementSym[nelem]) < 0 );
00467 }
00468
00469 if( sd.sdCase != SD_SINGLE_SIZE )
00470 {
00471 lgErr = lgErr || ( fprintf(fdes,"#\n# entire size distribution, amin=%.5f amax=%.5f micron\n",
00472 sd.lim[ipBLo],sd.lim[ipBHi]) < 0 );
00473 lgErr = lgErr || ( fprintf(fdes,"#\n%.6e # ratio a_max/a_min in each size bin\n",
00474 pow(sd.lim[ipBHi]/sd.lim[ipBLo],1./(double)sd.nPart) ) < 0 );
00475 lgErr = lgErr || ( fprintf(fdes,"#\n# size distribution function\n#\n") < 0 );
00476 lgErr = lgErr || ( fprintf(fdes,"%12ld # number of table entries\n#\n",NPTS_TABLE+1) < 0 );
00477 lgErr = lgErr || ( fprintf(fdes,"# ============================\n") < 0 );
00478 lgErr = lgErr || ( fprintf(fdes,"# size (micr) a^4*dN/da (cm^3/H)\n#\n") < 0 );
00479 for( i=0; i <= NPTS_TABLE; i++ )
00480 {
00481 double radius, a4dNda;
00482 radius = sd.lim[ipBLo]*exp((double)i/(double)NPTS_TABLE*log(sd.lim[ipBHi]/sd.lim[ipBLo]));
00483 radius = MAX2(MIN2(radius,sd.lim[ipBHi]),sd.lim[ipBLo]);
00484 a4dNda = POW4(radius)*size_distr(radius,&sd)/gd.norm*1.e-12/sd.unity;
00485 lgErr = lgErr || ( fprintf(fdes,"%.6e %.6e\n",radius,a4dNda) < 0 );
00486 }
00487 }
00488 else
00489 {
00490 lgErr = lgErr || ( fprintf(fdes,"#\n") < 0 );
00491 lgErr = lgErr || ( fprintf(fdes,"%.6e # a_max/a_min = 1 for single sized grain\n", 1. ) < 0 );
00492 lgErr = lgErr || ( fprintf(fdes,"%12ld # no size distribution table\n",0L) < 0 );
00493 }
00494
00495 lgErr = lgErr || ( fprintf(fdes,"#\n") < 0 );
00496 lgErr = lgErr || ( fprintf(fdes,"%12ld # rfield.nupper\n",rfield.nupper) < 0 );
00497 lgErr = lgErr || ( fprintf(fdes,"%12ld # number of size distr. bins\n#\n",sd.nPart) < 0 );
00498
00499 for( p=0; p < sd.nPart; p++ )
00500 {
00501 sd.cPart = p;
00502
00503 mie_auxiliary(&sd,"step");
00504
00505 if( sd.nPart > 1 )
00506 {
00507
00508
00509
00510
00511
00512 double frac = sd.unity_bin/sd.unity;
00513 volfrac = sd.vol*frac/volnorm;
00514 fprintf( ioQQQ, " Starting size bin %ld, amin=%.5f amax=%.5f micron\n",
00515 p+1,sd.clim[ipBLo],sd.clim[ipBHi] );
00516 lgErr = lgErr || ( fprintf(fdes,"# size bin %ld, amin=%.5f amax=%.5f micron\n",
00517 p+1,sd.clim[ipBLo],sd.clim[ipBHi]) < 0 );
00518 lgErr = lgErr || ( fprintf(fdes,"%.6e # average grain ",3.*sd.vol/sd.area) < 0 );
00519 lgErr = lgErr || ( fprintf(fdes,"radius <a^3>/<a^2>, this bin (cm)\n") < 0 );
00520 lgErr = lgErr || ( fprintf(fdes,"%.6e # average ",sd.area) < 0 );
00521 lgErr = lgErr || ( fprintf(fdes,"grain area <4pi*a^2>, this bin (cm^2)\n") < 0 );
00522 lgErr = lgErr || ( fprintf(fdes,"%.6e # average ",sd.vol) < 0 );
00523 lgErr = lgErr || ( fprintf(fdes,"grain volume <4/3pi*a^3>, this bin (cm^3)\n") < 0 );
00524 lgErr = lgErr || ( fprintf(fdes,"%.6e # total grain ",sd.radius*frac/gd.norm) < 0 );
00525 lgErr = lgErr || ( fprintf(fdes,"radius Int(a) per H, this bin (cm/H)\n") < 0 );
00526 lgErr = lgErr || ( fprintf(fdes,"%.6e # total grain area ",sd.area*frac/gd.norm) < 0 );
00527 lgErr = lgErr || ( fprintf(fdes,"Int(4pi*a^2) per H, this bin (cm^2/H)\n") < 0 );
00528 lgErr = lgErr || ( fprintf(fdes,"%.6e # total grain volume ",sd.vol*frac/gd.norm) < 0 );
00529 lgErr = lgErr || ( fprintf(fdes,"Int(4/3pi*a^3) per H, this bin (cm^3/H)\n#\n") < 0 );
00530 }
00531
00532 lgErrorOccurred = false;
00533
00534
00535 for( i=0; i < rfield.nupper; i++ )
00536 {
00537 wavlen = WAVNRYD/rfield.anu[i]*1.e4;
00538
00539 switch( gd.rfiType )
00540 {
00541 case RFI_TABLE:
00542 ErrorIndex[i] = 0;
00543 acs_abs[p][i] = 0.;
00544 acs_sct[p][i] = 0.;
00545 a1g[p][i] = 0.;
00546
00547 for( gd.cAxis=0; gd.cAxis < gd.nAxes; gd.cAxis++ )
00548 {
00549 mie_cs_size_distr(wavlen,&sd,&gd,mie_cs,&cs_abs,&cs_sct,&cosb,&Error);
00550 ErrorIndex[i] = MAX2(ErrorIndex[i],Error);
00551 acs_abs[p][i] += cs_abs*gd.wt[gd.cAxis];
00552 acs_sct[p][i] += cs_sct*gd.wt[gd.cAxis];
00553 a1g[p][i] += cs_sct*(1.-cosb)*gd.wt[gd.cAxis];
00554 }
00555
00556 if( ErrorIndex[i] > 0 )
00557 {
00558 ErrorIndex[i] = MIN2(ErrorIndex[i],2);
00559 lgErrorOccurred = true;
00560 }
00561
00562 switch( ErrorIndex[i] )
00563 {
00564
00565 case 2:
00566 acs_abs[p][i] = 0.;
00567 acs_sct[p][i] = 0.;
00568
00569
00570 case 1:
00571 a1g[p][i] = 0.;
00572 break;
00573
00574 case 0:
00575 a1g[p][i] /= acs_sct[p][i];
00576 break;
00577 default:
00578 fprintf( ioQQQ, " Insane value for ErrorIndex: %d\n", ErrorIndex[i] );
00579 ShowMe();
00580 puts( "[Stop in mie_write_opc]" );
00581 cdEXIT(EXIT_FAILURE);
00582 }
00583
00584
00585 if( ErrorIndex[i] < 2 )
00586 ASSERT( acs_abs[p][i] > 0. && acs_sct[p][i] > 0. );
00587 if( ErrorIndex[i] < 1 )
00588 ASSERT( a1g[p][i] > 0. );
00589 break;
00590 case OPC_TABLE:
00591
00592 gd.cAxis = 0;
00593 mie_cs_size_distr(wavlen,&sd,&gd,tbl_fun,&cs_abs,&cs_sct,&cosb,&Error);
00594 ErrorIndex[i] = MIN2(Error,2);
00595 lgErrorOccurred = lgErrorOccurred || ( Error > 0 );
00596 acs_abs[p][i] = cs_abs*gd.norm;
00597 acs_sct[p][i] = cs_sct*gd.norm;
00598 a1g[p][i] = 1.-cosb;
00599 break;
00600 case OPC_GREY:
00601 ErrorIndex[i] = 0;
00602 acs_abs[p][i] = 1.3121e-23*volfrac*gd.norm;
00603 acs_sct[p][i] = 2.6242e-23*volfrac*gd.norm;
00604 a1g[p][i] = 1.;
00605 break;
00606 case OPC_PAH1:
00607
00608 gd.cAxis = 0;
00609 mie_cs_size_distr(wavlen,&sd,&gd,pah1_fun,&cs_abs,&cs_sct,&cosb,&Error);
00610 ErrorIndex[i] = MIN2(Error,2);
00611 lgErrorOccurred = lgErrorOccurred || ( Error > 0 );
00612 acs_abs[p][i] = cs_abs;
00613 acs_sct[p][i] = cs_sct;
00614 a1g[p][i] = 1.-cosb;
00615 break;
00616 default:
00617 fprintf( ioQQQ, " Insanity in mie_write_opc\n" );
00618 ShowMe();
00619 puts( "[Stop in mie_write_opc]" );
00620 cdEXIT(EXIT_FAILURE);
00621 }
00622 }
00623
00624
00625 if( lgErrorOccurred )
00626 {
00627 strcpy(chString,"absorption cs");
00628 mie_repair(chString,rfield.nupper,2,0,rfield.anu,acs_abs[p],ErrorIndex,false,&lgWarning);
00629 strcpy(chString,"scattering cs");
00630 mie_repair(chString,rfield.nupper,2,1,rfield.anu,acs_sct[p],ErrorIndex,false,&lgWarning);
00631 strcpy(chString,"asymmetry parameter");
00632 mie_repair(chString,rfield.nupper,1,1,rfield.anu,a1g[p],ErrorIndex,true,&lgWarning);
00633 }
00634
00635 for( i=0; i < rfield.nupper; i++ )
00636 {
00637 acs_abs[p][i] /= gd.norm;
00638
00639
00640 acs_sct[p][i] /= gd.norm;
00641 }
00642 #if 0
00643 {
00644 FILE* ftmp;
00645 char name[20] = "asymxx";
00646 sprintf(&name[4],"%2.2li",p+1);
00647 ftmp = fopen(name,"w");
00648 for( i=0; i < rfield.nupper; i++ )
00649 {
00650 fprintf( ftmp, "%.6e %.6e\n", rfield.anu[i],a1g[p][i]);
00651 }
00652 fclose(ftmp);
00653 }
00654 #endif
00655
00656 mie_auxiliary(&sd,"cleanup");
00657 }
00658
00659 lgErr = lgErr || ( fprintf(fdes,"#\n") < 0 );
00660 lgErr = lgErr || ( fprintf(fdes,"# ===========================================\n") < 0 );
00661 lgErr = lgErr || ( fprintf(fdes,"# anu (Ryd) abs_cs_01 (cm^2/H) abs_cs_02.....\n#\n") < 0 );
00662
00663 for( i=0; i < rfield.nupper; i++ )
00664 {
00665 lgErr = lgErr || ( fprintf(fdes,"%.6e ",rfield.anu[i]) < 0 );
00666 for( p=0; p < sd.nPart; p++ )
00667 {
00668 lgErr = lgErr || ( fprintf(fdes,"%.6e ",acs_abs[p][i]) < 0 );
00669 }
00670 lgErr = lgErr || ( fprintf(fdes,"\n") < 0 );
00671 }
00672
00673 lgErr = lgErr || ( fprintf(fdes,"#\n") < 0 );
00674 lgErr = lgErr || ( fprintf(fdes,"# ===========================================\n") < 0 );
00675 lgErr = lgErr || ( fprintf(fdes,"# anu (Ryd) sct_cs_01 (cm^2/H) sct_cs_02.....\n#\n") < 0 );
00676
00677 for( i=0; i < rfield.nupper; i++ )
00678 {
00679 lgErr = lgErr || ( fprintf(fdes,"%.6e ",rfield.anu[i]) < 0 );
00680 for( p=0; p < sd.nPart; p++ )
00681 {
00682 lgErr = lgErr || ( fprintf(fdes,"%.6e ",acs_sct[p][i]) < 0 );
00683 }
00684 lgErr = lgErr || ( fprintf(fdes,"\n") < 0 );
00685 }
00686
00687 lgErr = lgErr || ( fprintf(fdes,"#\n") < 0 );
00688 lgErr = lgErr || ( fprintf(fdes,"# ===========================================\n") < 0 );
00689 lgErr = lgErr || ( fprintf(fdes,"# anu (Ryd) (1-g)_bin_01 (1-g)_bin_02.....\n#\n") < 0 );
00690
00691 for( i=0; i < rfield.nupper; i++ )
00692 {
00693 lgErr = lgErr || ( fprintf(fdes,"%.6e ",rfield.anu[i]) < 0 );
00694 for( p=0; p < sd.nPart; p++ )
00695 {
00696 lgErr = lgErr || ( fprintf(fdes,"%.6e ",a1g[p][i]) < 0 );
00697 }
00698 lgErr = lgErr || ( fprintf(fdes,"\n") < 0 );
00699 }
00700
00701 fprintf( ioQQQ, " Starting calculation of inverse attenuation length\n" );
00702 strcpy(chString,"inverse attenuation length");
00703 if( gd.rfiType != RFI_TABLE )
00704 {
00705
00706 ial_type icase = gv.which_ial[gd.matType];
00707 switch( icase )
00708 {
00709 case IAL_CAR:
00710 mie_read_rfi("graphite.rfi",&gd2);
00711 mie_calc_ial(&gd2,rfield.nupper,inv_att_len,chString,&lgWarning);
00712 break;
00713 case IAL_SIL:
00714 mie_read_rfi("silicate.rfi",&gd2);
00715 mie_calc_ial(&gd2,rfield.nupper,inv_att_len,chString,&lgWarning);
00716 break;
00717 default:
00718 fprintf( ioQQQ, " mie_write_opc detected unknown ial type: %d\n" , icase );
00719 puts( "[Stop in mie_write_opc]" );
00720 cdEXIT(EXIT_FAILURE);
00721 }
00722 }
00723 else
00724 {
00725 mie_calc_ial(&gd,rfield.nupper,inv_att_len,chString,&lgWarning);
00726 }
00727
00728 lgErr = lgErr || ( fprintf(fdes,"#\n") < 0 );
00729 lgErr = lgErr || ( fprintf(fdes,"# ===========================================\n") < 0 );
00730 lgErr = lgErr || ( fprintf(fdes,"# anu (Ryd) inverse attenuation length (cm^-1)\n#\n") < 0 );
00731
00732 for( i=0; i < rfield.nupper; i++ )
00733 {
00734 lgErr = lgErr || ( fprintf(fdes,"%.6e %.6e\n",rfield.anu[i],inv_att_len[i]) < 0 );
00735 }
00736
00737 fclose(fdes);
00738
00739 if( lgErr )
00740 {
00741 fprintf( ioQQQ, "\n Error writing file: %s\n", chFile );
00742 if( remove(chFile) == 0 )
00743 {
00744 fprintf( ioQQQ, " The file has been removed\n" );
00745 puts( "[Stop in mie_write_opc]" );
00746 cdEXIT(EXIT_FAILURE);
00747 }
00748 }
00749 else
00750 {
00751 fprintf( ioQQQ, "\n Opacity file %s written succesfully\n\n", chFile );
00752 if( lgWarning )
00753 {
00754 fprintf( ioQQQ, "\n !!! Warnings were detected !!!\n\n" );
00755 }
00756 }
00757
00758 for( p=0; p < sd.nPart; p++ )
00759 {
00760 free(a1g[p]);
00761 free(acs_sct[p]);
00762 free(acs_abs[p]);
00763 }
00764
00765 free(inv_att_len);
00766 free(a1g);
00767 free(acs_sct);
00768 free(acs_abs);
00769
00770
00771 if( sd.ln_a != NULL )
00772 free(sd.ln_a);
00773 if( sd.ln_a4dNda != NULL )
00774 free(sd.ln_a4dNda);
00775
00776
00777 for( j=0; j < gd2.nOpcCols; j++ )
00778 {
00779 if( gd2.opcData[j] != NULL )
00780 free(gd2.opcData[j]);
00781 }
00782 if( gd2.opcAnu != NULL )
00783 free(gd2.opcAnu);
00784 for( j=0; j < gd2.nAxes; j++ )
00785 {
00786 if( gd2.nr1[j] != NULL )
00787 free(gd2.nr1[j]);
00788 if( gd2.n[j] != NULL )
00789 free(gd2.n[j]);
00790 if( gd2.wavlen[j] != NULL )
00791 free(gd2.wavlen[j]);
00792 }
00793
00794 for( j=0; j < gd.nOpcCols; j++ )
00795 {
00796 if( gd.opcData[j] != NULL )
00797 free(gd.opcData[j]);
00798 }
00799 if( gd.opcAnu != NULL )
00800 free(gd.opcAnu);
00801 for( j=0; j < gd.nAxes; j++ )
00802 {
00803 if( gd.nr1[j] != NULL )
00804 free(gd.nr1[j]);
00805 if( gd.n[j] != NULL )
00806 free(gd.n[j]);
00807 if( gd.wavlen[j] != NULL )
00808 free(gd.wavlen[j]);
00809 }
00810
00811 free( ErrorIndex );
00812
00813 DEBUG_EXIT( "mie_write_opc()" );
00814 return;
00815 }
00816
00817
00818 static void mie_auxiliary( sd_data *sd,
00819 const char *auxCase)
00820 {
00821 double amin,
00822 amax,
00823 delta,
00824 oldvol,
00825 step;
00826
00827
00828 const double TOLER = 1.e-3;
00829
00830 DEBUG_ENTRY( "mie_auxiliary()" );
00831 if( strcmp(auxCase,"init") == 0 )
00832 {
00833 sd->xx = NULL;
00834 sd->aa = NULL;
00835 sd->rr = NULL;
00836 sd->ww = NULL;
00837
00838
00839
00840
00841 sd->nmul = 1;
00842
00843
00844 switch( sd->sdCase )
00845 {
00846 case SD_SINGLE_SIZE:
00847 sd->radius = sd->a[ipSize]*1.e-4;
00848 sd->area = 4.*PI*POW2(sd->a[ipSize])*1.e-8;
00849 sd->vol = 4./3.*PI*POW3(sd->a[ipSize])*1.e-12;
00850 break;
00851 case SD_POWERLAW:
00852 case SD_EXP_CUTOFF1:
00853 case SD_EXP_CUTOFF2:
00854 case SD_EXP_CUTOFF3:
00855 case SD_LOG_NORMAL:
00856 case SD_LIN_NORMAL:
00857 case SD_TABLE:
00858
00859
00860 amin = sd->lgLogScale ? log(sd->lim[ipBLo]) : sd->lim[ipBLo];
00861 amax = sd->lgLogScale ? log(sd->lim[ipBHi]) : sd->lim[ipBHi];
00862
00863 sd->clim[ipBLo] = sd->lim[ipBLo];
00864 sd->clim[ipBHi] = sd->lim[ipBHi];
00865
00866
00867 oldvol= 0.;
00868 do
00869 {
00870 sd->nmul *= 2;
00871 mie_integrate(sd,amin,amax,&sd->unity,true);
00872 delta = fabs(sd->vol-oldvol)/sd->vol;
00873 oldvol = sd->vol;
00874 } while( sd->nmul <= 1024 && delta > TOLER );
00875
00876 if( delta > TOLER )
00877 {
00878 fprintf( ioQQQ, " could not converge integration of size distribution\n" );
00879 puts( "[Stop in mie_auxiliary]" );
00880 cdEXIT(EXIT_FAILURE);
00881 }
00882
00883
00884
00885 sd->nmul /= 2;
00886 mie_integrate(sd,amin,amax,&sd->unity,true);
00887 break;
00888 case SD_ILLEGAL:
00889 default:
00890 fprintf( ioQQQ, " insane case for grain size distribution: %d\n" , sd->sdCase );
00891 ShowMe();
00892 puts( "[Stop in mie_auxiliary]" );
00893 cdEXIT(EXIT_FAILURE);
00894 }
00895 }
00896 else if( strcmp(auxCase,"step") == 0 )
00897 {
00898
00899 switch( sd->sdCase )
00900 {
00901 case SD_SINGLE_SIZE:
00902 break;
00903 case SD_POWERLAW:
00904 case SD_EXP_CUTOFF1:
00905 case SD_EXP_CUTOFF2:
00906 case SD_EXP_CUTOFF3:
00907 case SD_LOG_NORMAL:
00908 case SD_LIN_NORMAL:
00909 case SD_TABLE:
00910 amin = sd->lgLogScale ? log(sd->lim[ipBLo]) : sd->lim[ipBLo];
00911 amax = sd->lgLogScale ? log(sd->lim[ipBHi]) : sd->lim[ipBHi];
00912 step = (amax - amin)/(double)sd->nPart;
00913 amin = amin + (double)sd->cPart*step;
00914 amax = MIN2(amax,amin + step);
00915
00916 sd->clim[ipBLo] = sd->lgLogScale ? exp(amin) : amin;
00917 sd->clim[ipBHi] = sd->lgLogScale ? exp(amax) : amax;
00918
00919 mie_integrate(sd,amin,amax,&sd->unity_bin,false);
00920
00921 break;
00922 case SD_ILLEGAL:
00923 default:
00924 fprintf( ioQQQ, " insane case for grain size distribution: %d\n" , sd->sdCase );
00925 ShowMe();
00926 puts( "[Stop in mie_auxiliary]" );
00927 cdEXIT(EXIT_FAILURE);
00928 }
00929 }
00930 else if( strcmp(auxCase,"cleanup") == 0 )
00931 {
00932 if( sd->xx != NULL )
00933 free(sd->xx);
00934 sd->xx = NULL;
00935 if( sd->aa != NULL )
00936 free(sd->aa);
00937 sd->aa = NULL;
00938 if( sd->rr != NULL )
00939 free(sd->rr);
00940 sd->rr = NULL;
00941 if( sd->ww != NULL )
00942 free(sd->ww);
00943 sd->ww = NULL;
00944 }
00945 else
00946 {
00947 fprintf( ioQQQ, " mie_auxiliary called with insane argument: %s\n", auxCase );
00948 ShowMe();
00949 puts( "[Stop in mie_auxiliary]" );
00950 cdEXIT(EXIT_FAILURE);
00951 }
00952
00953 DEBUG_EXIT( "mie_auxiliary()" );
00954 return;
00955 }
00956
00957 static void mie_integrate( sd_data *sd,
00958 double amin,
00959 double amax,
00960 double *normalization,
00961 bool lgFreeMem)
00962 {
00963 long int j;
00964 double unity;
00965
00966 DEBUG_ENTRY( "mie_integrate()" );
00967
00968
00969
00970 sd->nn = sd->nmul*((long)(2.*log(sd->clim[ipBHi]/sd->clim[ipBLo])) + 1);
00971 sd->nn = MIN2(MAX2(sd->nn,2*sd->nmul),4096);
00972 sd->xx = (double *)MALLOC(sizeof(double)*(unsigned)sd->nn);
00973 sd->aa = (double *)MALLOC(sizeof(double)*(unsigned)sd->nn);
00974 sd->rr = (double *)MALLOC(sizeof(double)*(unsigned)sd->nn);
00975 sd->ww = (double *)MALLOC(sizeof(double)*(unsigned)sd->nn);
00976 if( sd->xx == NULL || sd->aa == NULL || sd->rr == NULL || sd->ww == NULL )
00977 {
00978 fprintf( ioQQQ, " Could not MALLOC arrays\n" );
00979 puts( "[Stop in mie_integrate]" );
00980 cdEXIT(EXIT_FAILURE);
00981 }
00982 gauss_legendre(sd->nn,sd->xx,sd->aa);
00983 gauss_init(sd->nn,amin,amax,sd->xx,sd->aa,sd->rr,sd->ww);
00984
00985
00986 unity = 0.;
00987 sd->radius = 0.;
00988 sd->area = 0.;
00989 sd->vol = 0.;
00990
00991 for( j=0; j < sd->nn; j++ )
00992 {
00993 double weight;
00994
00995
00996 if( sd->lgLogScale )
00997 {
00998 sd->rr[j] = exp(sd->rr[j]);
00999 sd->ww[j] *= sd->rr[j];
01000 }
01001 weight = sd->ww[j]*size_distr(sd->rr[j],sd);
01002 unity += weight;
01003 sd->radius += weight*sd->rr[j];
01004 sd->area += weight*POW2(sd->rr[j]);
01005 sd->vol += weight*POW3(sd->rr[j]);
01006 }
01007 *normalization = unity;
01008 sd->radius *= 1.e-4/unity;
01009 sd->area *= 4.*PI*1.e-8/unity;
01010 sd->vol *= 4./3.*PI*1.e-12/unity;
01011
01012 if( lgFreeMem )
01013 {
01014 if( sd->xx != NULL )
01015 free(sd->xx);
01016 sd->xx = NULL;
01017 if( sd->aa != NULL )
01018 free(sd->aa);
01019 sd->aa = NULL;
01020 if( sd->rr != NULL )
01021 free(sd->rr);
01022 sd->rr = NULL;
01023 if( sd->ww != NULL )
01024 free(sd->ww);
01025 sd->ww = NULL;
01026 }
01027
01028 DEBUG_EXIT( "mie_integrate()" );
01029 return;
01030 }
01031
01032
01033 void mie_read_opc(const char *chFile,
01034 GrainPar gp)
01035 {
01036 int res,
01037 lgDefaultQHeat;
01038 long int dl,
01039 help,
01040 i,
01041 nelem,
01042 j,
01043 magic,
01044 nbin,
01045 nd,
01046 nd2,
01047 nup;
01048 float RefAbund[LIMELM],
01049 VolTotal;
01050 double anu;
01051 double RadiusRatio;
01052 char chLine[FILENAME_PATH_LENGTH_2],
01053 *str,
01054 chString[FILENAME_PATH_LENGTH_2];
01055 FILE *io2;
01056
01057
01058
01059 const double RATIO_MAX = pow(100.,1./3.);
01060
01061 DEBUG_ENTRY( "mie_read_opc()" );
01062
01063
01064 if( lgDataPathSet )
01065 {
01066
01067 strcpy( chString , chDataPath );
01068 strcat( chString, chFile );
01069 }
01070 else
01071 {
01072
01073 strcpy( chString , chFile );
01074 }
01075
01076 if( (io2 = fopen(chString,"r")) == NULL )
01077 {
01078
01079 strcpy( chString , chFile );
01080 if( (io2 = fopen(chString,"r")) == NULL )
01081 {
01082 fprintf( ioQQQ, " Could not open %s for reading\n",chString);
01083
01084 path_not_set();
01085 puts( "[Stop in mie_read_opc]" );
01086 cdEXIT(EXIT_FAILURE);
01087 }
01088 }
01089
01090
01091 strcpy( chLine, " " );
01092 if( strlen(chFile) <= 40 )
01093 {
01094 strncpy( &chLine[0], chFile, strlen(chFile) );
01095 }
01096 else
01097 {
01098 strncpy( &chLine[0], chFile, 37 );
01099 strncpy( &chLine[37], "...", 3 );
01100 }
01101 if( called.lgTalk )
01102 fprintf( ioQQQ, " * >>>> mie_read_opc reading file -- %40s <<<< *\n", chLine );
01103
01104
01105 for( i=0; i < gv.ReadPtr; ++i )
01106 {
01107 if( strcmp(chFile,gv.ReadRecord[i]) == 0 )
01108 {
01109 fprintf( ioQQQ, " File %s has already been read before, was this intended ?\n", chFile );
01110 break;
01111 }
01112 }
01113
01114 if( gv.ReadPtr < MAX_READ_RECORDS )
01115 {
01116 strcpy(gv.ReadRecord[gv.ReadPtr],chFile);
01117 ++gv.ReadPtr;
01118 }
01119
01120
01121 nd = NewGrainBin();
01122
01123 dl = 0;
01124
01125
01126 mie_next_data(chFile,io2,chLine,&dl);
01127 mie_read_long(chFile,chLine,&magic,true,dl);
01128 if( magic != MAGIC_OPC )
01129 {
01130 fprintf( ioQQQ, " Opacity file %s has obsolete magic number\n",chFile );
01131 fprintf( ioQQQ, " I found magic number %ld, but expected %ld on line #%ld\n",
01132 magic,MAGIC_OPC,dl );
01133 fprintf( ioQQQ, " Please recompile this file\n" );
01134 puts( "[Stop in mie_read_opc]" );
01135 cdEXIT(EXIT_FAILURE);
01136 }
01137
01138
01139 mie_next_data(chFile,io2,chLine,&dl);
01140 mie_read_long(chFile,chLine,&magic,true,dl);
01141
01142 mie_next_data(chFile,io2,chLine,&dl);
01143 mie_read_long(chFile,chLine,&magic,true,dl);
01144
01145
01146
01147
01148
01149
01150 gv.bin[nd]->dstfactor = (float)gp.dep;
01151
01152
01153 mie_next_data(chFile,io2,chLine,&dl);
01154 strncpy(gv.bin[nd]->chDstLab,chLine,(size_t)LABELSIZE);
01155 gv.bin[nd]->chDstLab[LABELSIZE] = '\0';
01156
01157
01158 mie_next_data(chFile,io2,chLine,&dl);
01159 mie_read_float(chFile,chLine,&gv.bin[nd]->dustp[0],true,dl);
01160
01161
01162 mie_next_data(chFile,io2,chLine,&dl);
01163 mie_read_float(chFile,chLine,&gv.bin[nd]->dustp[1],true,dl);
01164
01165
01166 mie_next_data(chFile,io2,chLine,&dl);
01167 mie_read_float(chFile,chLine,&gv.bin[nd]->atomWeight,true,dl);
01168
01169
01170 mie_next_data(chFile,io2,chLine,&dl);
01171 mie_read_float(chFile,chLine,&gv.bin[nd]->dustp[2],true,dl);
01172
01173
01174 mie_next_data(chFile,io2,chLine,&dl);
01175 mie_read_float(chFile,chLine,&gv.bin[nd]->dustp[3],true,dl);
01176
01177
01178 mie_next_data(chFile,io2,chLine,&dl);
01179 mie_read_float(chFile,chLine,&gv.bin[nd]->AvRadius,true,dl);
01180
01181
01182 mie_next_data(chFile,io2,chLine,&dl);
01183 mie_read_float(chFile,chLine,&gv.bin[nd]->AvArea,true,dl);
01184
01185
01186 mie_next_data(chFile,io2,chLine,&dl);
01187 mie_read_float(chFile,chLine,&gv.bin[nd]->AvVol,true,dl);
01188
01189
01190 mie_next_data(chFile,io2,chLine,&dl);
01191 mie_read_float(chFile,chLine,&gv.bin[nd]->IntRadius,true,dl);
01192
01193
01194 mie_next_data(chFile,io2,chLine,&dl);
01195 mie_read_float(chFile,chLine,&gv.bin[nd]->IntArea,true,dl);
01196
01197
01198 mie_next_data(chFile,io2,chLine,&dl);
01199 mie_read_float(chFile,chLine,&gv.bin[nd]->IntVol,true,dl);
01200
01201
01202 mie_next_data(chFile,io2,chLine,&dl);
01203 mie_read_float(chFile,chLine,&gv.bin[nd]->DustWorkFcn,true,dl);
01204
01205
01206 mie_next_data(chFile,io2,chLine,&dl);
01207 mie_read_float(chFile,chLine,&gv.bin[nd]->BandGap,false,dl);
01208
01209
01210 mie_next_data(chFile,io2,chLine,&dl);
01211 mie_read_float(chFile,chLine,&gv.bin[nd]->ThermEff,true,dl);
01212
01213
01214 mie_next_data(chFile,io2,chLine,&dl);
01215 mie_read_float(chFile,chLine,&gv.bin[nd]->Tsublimat,true,dl);
01216
01217
01218 mie_next_data(chFile,io2,chLine,&dl);
01219 mie_read_long(chFile,chLine,&help,true,dl);
01220 gv.bin[nd]->matType = (mat_type)help;
01221
01222 for( nelem=0; nelem < LIMELM; nelem++ )
01223 {
01224 mie_next_data(chFile,io2,chLine,&dl);
01225 mie_read_float(chFile,chLine,&RefAbund[nelem],false,dl);
01226
01227
01228 gv.bin[nd]->AccomCoef[nelem] = 2.*gv.bin[nd]->atomWeight*dense.AtomicWeight[nelem]/
01229 POW2(gv.bin[nd]->atomWeight+dense.AtomicWeight[nelem]);
01230 }
01231
01232
01233 mie_next_data(chFile,io2,chLine,&dl);
01234 mie_read_double(chFile,chLine,&RadiusRatio,true,dl);
01235
01236 gv.bin[nd]->lgDustFunc = gp.lgAbunVsDepth;
01237
01238 lgDefaultQHeat = ( RadiusRatio < RATIO_MAX && !gp.lgGreyGrain );
01239 gv.bin[nd]->lgQHeat = ( gp.lgForbidQHeating ) ? false : ( gp.lgRequestQHeating || lgDefaultQHeat );
01240 gv.bin[nd]->cnv_H_pGR = gv.bin[nd]->AvVol/gv.bin[nd]->IntVol;
01241 gv.bin[nd]->cnv_GR_pH = 1./gv.bin[nd]->cnv_H_pGR;
01242
01243
01244 gv.bin[nd]->Capacity = PI4*ELECTRIC_CONST*gv.bin[nd]->IntRadius/100.*gv.bin[nd]->cnv_H_pGR;
01245
01246
01247 mie_next_data(chFile,io2,chLine,&dl);
01248 mie_read_long(chFile,chLine,&nup,false,dl);
01249 for( i=0; i < nup; i++ )
01250 mie_next_data(chFile,io2,chLine,&dl);
01251
01252
01253 mie_next_data(chFile,io2,chLine,&dl);
01254 mie_read_long(chFile,chLine,&nup,true,dl);
01255
01256 gv.bin[nd]->NFPCheck = nup;
01257
01258
01259 mie_next_data(chFile,io2,chLine,&dl);
01260 mie_read_long(chFile,chLine,&nbin,true,dl);
01261
01262 if( nbin == 1 )
01263 {
01264
01265 ASSERT( gv.bin[nd]->dstab1 == NULL );
01266 gv.bin[nd]->dstab1 = (double*)MALLOC(sizeof(double)*(unsigned)nup);
01267 ASSERT( gv.bin[nd]->pure_sc1 == NULL );
01268 gv.bin[nd]->pure_sc1 = (double*)MALLOC(sizeof(double)*(unsigned)nup);
01269 ASSERT( gv.bin[nd]->asym == NULL );
01270 gv.bin[nd]->asym = (double*)MALLOC(sizeof(double)*(unsigned)nup);
01271 ASSERT( gv.bin[nd]->inv_att_len == NULL );
01272 gv.bin[nd]->inv_att_len = (float*)MALLOC(sizeof(float)*(unsigned)nup);
01273
01274 gv.bin[nd]->dustp[4] = 1.;
01275 for( nelem=0; nelem < LIMELM; nelem++ )
01276 {
01277 gv.bin[nd]->elmAbund[nelem] = RefAbund[nelem];
01278 }
01279 }
01280 else if( nbin > 1 )
01281 {
01282
01283 VolTotal = gv.bin[nd]->IntVol;
01284
01285 for( i=0; i < nbin; i++ )
01286 {
01287
01288 nd2 = ( i == 0 ) ? nd : NewGrainBin();
01289
01290
01291 ASSERT( gv.bin[nd2]->dstab1 == NULL );
01292 gv.bin[nd2]->dstab1 = (double*)MALLOC(sizeof(double)*(unsigned)nup);
01293 ASSERT( gv.bin[nd2]->pure_sc1 == NULL );
01294 gv.bin[nd2]->pure_sc1 = (double*)MALLOC(sizeof(double)*(unsigned)nup);
01295 ASSERT( gv.bin[nd2]->asym == NULL );
01296 gv.bin[nd2]->asym = (double*)MALLOC(sizeof(double)*(unsigned)nup);
01297 ASSERT( gv.bin[nd2]->inv_att_len == NULL );
01298 gv.bin[nd2]->inv_att_len = (float*)MALLOC(sizeof(float)*(unsigned)nup);
01299
01300
01301 mie_next_data(chFile,io2,chLine,&dl);
01302 mie_read_float(chFile,chLine,&gv.bin[nd2]->AvRadius,true,dl);
01303
01304
01305 mie_next_data(chFile,io2,chLine,&dl);
01306 mie_read_float(chFile,chLine,&gv.bin[nd2]->AvArea,true,dl);
01307
01308
01309 mie_next_data(chFile,io2,chLine,&dl);
01310 mie_read_float(chFile,chLine,&gv.bin[nd2]->AvVol,true,dl);
01311
01312
01313 mie_next_data(chFile,io2,chLine,&dl);
01314 mie_read_float(chFile,chLine,&gv.bin[nd2]->IntRadius,true,dl);
01315
01316
01317 mie_next_data(chFile,io2,chLine,&dl);
01318 mie_read_float(chFile,chLine,&gv.bin[nd2]->IntArea,true,dl);
01319
01320
01321 mie_next_data(chFile,io2,chLine,&dl);
01322 mie_read_float(chFile,chLine,&gv.bin[nd2]->IntVol,true,dl);
01323
01324 gv.bin[nd2]->cnv_H_pGR = gv.bin[nd2]->AvVol/gv.bin[nd2]->IntVol;
01325 gv.bin[nd2]->cnv_GR_pH = 1./gv.bin[nd2]->cnv_H_pGR;
01326
01327
01328
01329
01330 gv.bin[nd2]->Capacity =
01331 PI4*ELECTRIC_CONST*gv.bin[nd2]->IntRadius/100.*gv.bin[nd2]->cnv_H_pGR;
01332
01333
01334
01335
01336 gv.bin[nd2]->dustp[4] = gv.bin[nd2]->IntVol/VolTotal;
01337 for( nelem=0; nelem < LIMELM; nelem++ )
01338 {
01339 gv.bin[nd2]->elmAbund[nelem] = RefAbund[nelem]*gv.bin[nd2]->dustp[4];
01340 }
01341
01342 if( i > 0 )
01343 {
01344 gv.bin[nd2]->dstfactor = gv.bin[nd]->dstfactor;
01345 strcpy(gv.bin[nd2]->chDstLab,gv.bin[nd]->chDstLab);
01346 gv.bin[nd2]->atomWeight = gv.bin[nd]->atomWeight;
01347 for( j=0; j < 4; j++ )
01348 {
01349 gv.bin[nd2]->dustp[j] = gv.bin[nd]->dustp[j];
01350 }
01351 gv.bin[nd2]->DustWorkFcn = gv.bin[nd]->DustWorkFcn;
01352 gv.bin[nd2]->BandGap = gv.bin[nd]->BandGap;
01353 gv.bin[nd2]->ThermEff = gv.bin[nd]->ThermEff;
01354 gv.bin[nd2]->Tsublimat = gv.bin[nd]->Tsublimat;
01355 gv.bin[nd2]->matType = gv.bin[nd]->matType;
01356 gv.bin[nd2]->lgDustFunc = gv.bin[nd]->lgDustFunc;
01357 gv.bin[nd2]->lgQHeat = gv.bin[nd]->lgQHeat;
01358 gv.bin[nd2]->NFPCheck = gv.bin[nd]->NFPCheck;
01359 for( nelem=0; nelem < LIMELM; nelem++ )
01360 {
01361 gv.bin[nd2]->AccomCoef[nelem] = gv.bin[nd]->AccomCoef[nelem];
01362 }
01363 }
01364 }
01365 for( i=0; i < nbin; i++ )
01366 {
01367 nd2 = nd + i;
01368
01369 str = strstr(gv.bin[nd2]->chDstLab,"xx");
01370 if( str != NULL )
01371 sprintf(str,"%02ld",i+1);
01372 }
01373 }
01374
01375
01376 for( i=0; i < 5; i++ )
01377 mie_next_line(chFile,io2,chLine,&dl);
01378
01379
01380 for( i=0; i < nup; i++ )
01381 {
01382
01383 if( (res = fscanf(io2,"%le",&anu)) != 1 )
01384 {
01385 fprintf( ioQQQ, " Read failed on %s\n",chFile );
01386 if( res == EOF )
01387 fprintf( ioQQQ, " EOF reached prematurely\n" );
01388 puts( "[Stop in mie_read_opc]" );
01389 cdEXIT(EXIT_FAILURE);
01390 }
01391 for( j=0; j < nbin; j++ )
01392 {
01393 nd2 = nd + j;
01394 if( (res = fscanf(io2,"%le",&gv.bin[nd2]->dstab1[i])) != 1 )
01395 {
01396 fprintf( ioQQQ, " Read failed on %s\n",chFile );
01397 if( res == EOF )
01398 fprintf( ioQQQ, " EOF reached prematurely\n" );
01399 puts( "[Stop in mie_read_opc]" );
01400 cdEXIT(EXIT_FAILURE);
01401 }
01402 ASSERT( gv.bin[nd2]->dstab1[i] > 0. );
01403 }
01404 }
01405
01406
01407 for( i=0; i < 5; i++ )
01408 mie_next_line(chFile,io2,chLine,&dl);
01409
01410
01411 for( i=0; i < nup; i++ )
01412 {
01413 if( (res = fscanf(io2,"%le",&anu)) != 1 )
01414 {
01415 fprintf( ioQQQ, " Read failed on %s\n",chFile );
01416 if( res == EOF )
01417 fprintf( ioQQQ, " EOF reached prematurely\n" );
01418 puts( "[Stop in mie_read_opc]" );
01419 cdEXIT(EXIT_FAILURE);
01420 }
01421 for( j=0; j < nbin; j++ )
01422 {
01423 nd2 = nd + j;
01424 if( (res = fscanf(io2,"%le",&gv.bin[nd2]->pure_sc1[i])) != 1 )
01425 {
01426 fprintf( ioQQQ, " Read failed on %s\n",chFile );
01427 if( res == EOF )
01428 fprintf( ioQQQ, " EOF reached prematurely\n" );
01429 puts( "[Stop in mie_read_opc]" );
01430 cdEXIT(EXIT_FAILURE);
01431 }
01432 ASSERT( gv.bin[nd2]->pure_sc1[i] > 0. );
01433 }
01434 }
01435
01436
01437
01438
01439 for( i=0; i < 5; i++ )
01440 mie_next_line(chFile,io2,chLine,&dl);
01441
01442
01443 for( i=0; i < nup; i++ )
01444 {
01445 if( (res = fscanf(io2,"%le",&anu)) != 1 )
01446 {
01447 fprintf( ioQQQ, " Read failed on %s\n",chFile );
01448 if( res == EOF )
01449 fprintf( ioQQQ, " EOF reached prematurely\n" );
01450 puts( "[Stop in mie_read_opc]" );
01451 cdEXIT(EXIT_FAILURE);
01452 }
01453 for( j=0; j < nbin; j++ )
01454 {
01455 nd2 = nd + j;
01456 if( (res = fscanf(io2,"%le",&gv.bin[nd2]->asym[i])) != 1 )
01457 {
01458 fprintf( ioQQQ, " Read failed on %s\n",chFile );
01459 if( res == EOF )
01460 fprintf( ioQQQ, " EOF reached prematurely\n" );
01461 puts( "[Stop in mie_read_opc]" );
01462 cdEXIT(EXIT_FAILURE);
01463 }
01464 ASSERT( gv.bin[nd2]->asym[i] > 0. );
01465 }
01466 }
01467
01468
01469 for( i=0; i < 5; i++ )
01470 mie_next_line(chFile,io2,chLine,&dl);
01471
01472
01473 for( i=0; i < nup; i++ )
01474 {
01475 if( (res = fscanf(io2,"%le %e",&anu,&gv.bin[nd]->inv_att_len[i])) != 2 )
01476 {
01477 fprintf( ioQQQ, " Read failed on %s\n",chFile );
01478 if( res == EOF )
01479 fprintf( ioQQQ, " EOF reached prematurely\n" );
01480 puts( "[Stop in mie_read_opc]" );
01481 cdEXIT(EXIT_FAILURE);
01482 }
01483 ASSERT( gv.bin[nd]->inv_att_len[i] > 0. );
01484
01485 gv.bin[nd]->EnergyCheck = (float)anu;
01486
01487 for( j=1; j < nbin; j++ )
01488 {
01489 nd2 = nd + j;
01490 gv.bin[nd2]->inv_att_len[i] = gv.bin[nd]->inv_att_len[i];
01491 gv.bin[nd2]->EnergyCheck = gv.bin[nd]->EnergyCheck;
01492 }
01493 }
01494
01495 fclose(io2);
01496
01497 DEBUG_EXIT( "mie_read_opc()" );
01498 return;
01499 }
01500
01501
01502
01503
01504 static void mie_cs_size_distr(double wavlen,
01505 sd_data *sd,
01506 grain_data *gd,
01507 void(*cs_fun)(double,sd_data*,grain_data*,
01508 double*,double*,
01509 double*,int*),
01510 double *cs_abs,
01511 double *cs_sct,
01512 double *cosb,
01513 int *error)
01514 {
01515 bool lgADLused;
01516 long int i;
01517 double absval,
01518 g,
01519 sctval,
01520 weight;
01521
01522 DEBUG_ENTRY( "mie_cs_size_distr()" );
01523
01524
01525 ASSERT( wavlen > 0. );
01526 ASSERT( gd->cAxis >= 0 && gd->cAxis < gd->nAxes && gd->cAxis < NAX );
01527
01528 switch( sd->sdCase )
01529 {
01530 case SD_SINGLE_SIZE:
01531
01532 ASSERT( sd->a[ipSize] > 0. );
01533 sd->cSize = sd->a[ipSize];
01534 (*cs_fun)(wavlen,sd,gd,cs_abs,cs_sct,cosb,error);
01535 break;
01536 case SD_POWERLAW:
01537
01538 case SD_EXP_CUTOFF1:
01539 case SD_EXP_CUTOFF2:
01540 case SD_EXP_CUTOFF3:
01541
01542 case SD_LOG_NORMAL:
01543
01544 case SD_LIN_NORMAL:
01545
01546 case SD_TABLE:
01547
01548 ASSERT( sd->lim[ipBLo] > 0. && sd->lim[ipBHi] > 0. && sd->lim[ipBHi] > sd->lim[ipBLo] );
01549 lgADLused = false;
01550 *cs_abs = 0.;
01551 *cs_sct = 0.;
01552 *cosb = 0.;
01553 for( i=0; i < sd->nn; i++ )
01554 {
01555 sd->cSize = sd->rr[i];
01556 (*cs_fun)(wavlen,sd,gd,&absval,&sctval,&g,error);
01557 if( *error >= 2 )
01558 {
01559
01560 *cs_abs = -1.;
01561 *cs_sct = -1.;
01562 *cosb = -2.;
01563
01564 DEBUG_EXIT( "mie_cs_size_distr()" );
01565 return;
01566 }
01567 else if( *error == 1 )
01568 {
01569
01570 lgADLused = true;
01571 }
01572 weight = sd->ww[i]*size_distr(sd->rr[i],sd);
01573 *cs_abs += weight*absval;
01574 *cs_sct += weight*sctval;
01575 *cosb += weight*sctval*g;
01576 }
01577 if( lgADLused )
01578 {
01579 *error = 1;
01580 *cosb = -2.;
01581 }
01582 else
01583 {
01584 *error = 0;
01585 *cosb /= *cs_sct;
01586 }
01587 *cs_abs /= sd->unity;
01588 *cs_sct /= sd->unity;
01589 break;
01590 case SD_ILLEGAL:
01591 default:
01592 fprintf( ioQQQ, " insane case for grain size distribution: %d\n" , sd->sdCase );
01593 ShowMe();
01594 puts( "[Stop in mie_cs_size_distr]" );
01595 cdEXIT(EXIT_FAILURE);
01596 }
01597
01598 if( *error < 2 )
01599 ASSERT( *cs_abs > 0. && *cs_sct > 0. );
01600 if( *error < 1 )
01601 ASSERT( fabs(*cosb) <= 1.+10.*DBL_EPSILON );
01602
01603 DEBUG_EXIT( "mie_cs_size_distr()" );
01604 return;
01605 }
01606
01607
01608
01609
01610 static void mie_cs(double wavlen,
01611 sd_data *sd,
01612 grain_data *gd,
01613 double *cs_abs,
01614 double *cs_sct,
01615 double *cosb,
01616 int *error)
01617 {
01618 bool lgOutOfBounds;
01619 long int iflag,
01620 ind;
01621 double area,
01622 aqabs,
01623 aqext,
01624 aqphas,
01625 beta,
01626 ctbrqs = -DBL_MAX,
01627 delta,
01628 frac,
01629 nim,
01630 nre,
01631 nr1,
01632 qback,
01633 qext = -DBL_MAX,
01634 qphase,
01635 qscatt = -DBL_MAX,
01636 x,
01637 xistar;
01638
01639 DEBUG_ENTRY( "mie_cs()" );
01640
01641
01642 ASSERT( wavlen > 0. );
01643 ASSERT( sd->cSize > 0. );
01644 ASSERT( gd->wavlen[gd->cAxis] != NULL && gd->ndata[gd->cAxis] > 1 );
01645
01646
01647
01648
01649 find_arr(wavlen,gd->wavlen[gd->cAxis],gd->ndata[gd->cAxis],&ind,&lgOutOfBounds);
01650
01651 if( lgOutOfBounds )
01652 {
01653 *error = 3;
01654 *cs_abs = -1.;
01655 *cs_sct = -1.;
01656 *cosb = -2.;
01657
01658 DEBUG_EXIT( "mie_cs()" );
01659 return;
01660 }
01661
01662 frac = (wavlen-gd->wavlen[gd->cAxis][ind])/(gd->wavlen[gd->cAxis][ind+1]-gd->wavlen[gd->cAxis][ind]);
01663 ASSERT( frac > 0.-10.*DBL_EPSILON && frac < 1.+10.*DBL_EPSILON );
01664 nre = (1.-frac)*gd->n[gd->cAxis][ind].real() + frac*gd->n[gd->cAxis][ind+1].real();
01665 ASSERT( nre > 0. );
01666 nim = (1.-frac)*gd->n[gd->cAxis][ind].imag() + frac*gd->n[gd->cAxis][ind+1].imag();
01667 ASSERT( nim > 0. );
01668 nr1 = (1.-frac)*gd->nr1[gd->cAxis][ind] + frac*gd->nr1[gd->cAxis][ind+1];
01669 ASSERT( fabs(nre-1.-nr1) < 10.*MAX2(nre,1.)*DBL_EPSILON );
01670
01671
01672 area = PI*POW2(sd->cSize)*1.e-8;
01673
01674 x = sd->cSize/wavlen*2.*PI;
01675
01676
01677
01678
01679 sinpar(nre,nim,x,&qext,&qphase,&qscatt,&ctbrqs,&qback,&iflag);
01680
01681
01682
01683
01684 if( iflag == 0 )
01685 {
01686 *error = 0;
01687 *cs_abs = area*(qext - qscatt);
01688 *cs_sct = area*qscatt;
01689 *cosb = ctbrqs/qscatt;
01690 }
01691 else
01692 {
01693
01694 if( x >= 100. && sqrt(nr1*nr1+nim*nim) <= 0.001 )
01695 {
01696 delta = -nr1;
01697 beta = nim;
01698
01699 anomal(x,&aqext,&aqabs,&aqphas,&xistar,delta,beta);
01700
01701
01702 *error = 1;
01703 *cs_abs = area*aqabs;
01704 *cs_sct = area*(aqext - aqabs);
01705 *cosb = -2.;
01706 }
01707
01708 else
01709 {
01710 *error = 2;
01711 *cs_abs = -1.;
01712 *cs_sct = -1.;
01713 *cosb = -2.;
01714 }
01715 }
01716 if( *error < 2 )
01717 {
01718 if( *cs_abs <= 0. || *cs_sct <= 0. )
01719 {
01720 fprintf( ioQQQ, " illegal opacity found: wavl=%.4e micron," , wavlen );
01721 fprintf( ioQQQ, " abs_cs=%.2e, sct_cs=%.2e\n" , *cs_abs , *cs_sct );
01722 fprintf( ioQQQ, " please check refractive index file...\n" );
01723 puts( "[Stop in mie_cs]" );
01724 cdEXIT(EXIT_FAILURE);
01725 }
01726 }
01727 if( *error < 1 )
01728 {
01729 if( fabs(*cosb) > 1.+10.*DBL_EPSILON )
01730 {
01731 fprintf( ioQQQ, " illegal asymmetry parameter found: wavl=%.4e micron," , wavlen );
01732 fprintf( ioQQQ, " cosb=%.2e\n" , *cosb );
01733 fprintf( ioQQQ, " please check refractive index file...\n" );
01734 puts( "[Stop in mie_cs]" );
01735 cdEXIT(EXIT_FAILURE);
01736 }
01737 }
01738
01739
01740 DEBUG_EXIT( "mie_cs()" );
01741 return;
01742 }
01743
01744
01745
01746
01747
01748
01749
01750 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};
01751 static const double pah1_wlBand[7]={3.3, 6.18, 7.7, 8.6, 11.3, 12.0, 13.3};
01752 static const double pah1_width[7]={0.024, 0.102, 0.24, 0.168, 0.086, 0.174, 0.174};
01753
01754 static void pah1_fun(double wavl,
01755 sd_data *sd,
01756 grain_data *gd,
01757 double *cs_abs,
01758 double *cs_sct,
01759 double *cosb,
01760 int *error)
01761 {
01762 long int j;
01763 double cval1,
01764 pah1_fun_v,
01765 term,
01766 term1,
01767 term2,
01768 term3,
01769 x;
01770
01771 const double p1 = 4.0e-18;
01772 const double p2 = 1.1e-18;
01773 const double p3 = 3.3e-21;
01774 const double p4 = 6.0e-21;
01775 const double p5 = 2.4e-21;
01776 const double wl1a = 5.0;
01777 const double wl1b = 7.0;
01778 const double wl1c = 9.0;
01779 const double wl1d = 9.5;
01780 const double wl2a = 11.0;
01781 const double delwl2 = 0.3;
01782
01783 const double wl2b = wl2a + delwl2;
01784 const double wl2c = 15.0;
01785 const double wl3a = 3.2;
01786 const double wl3b = 3.57;
01787 const double wl3m = (wl3a+wl3b)/2.;
01788 const double wl3sig = 0.1476;
01789 const double x1 = 4.0;
01790 const double x2 = 5.9;
01791 const double x2a = RYD_INF/1.e4;
01792 const double x3 = 0.1;
01793
01794
01795 double vol = 4.*PI/3.*POW3(sd->cSize)*1.e-12;
01796
01797 double xnc = floor(vol*gd->rho/(ATOMIC_MASS_UNIT*dense.AtomicWeight[ipCARBON]));
01798
01799
01800 double xnh = floor(sqrt(6.*xnc));
01801
01802 double xnhoc = xnh/xnc;
01803
01804 double ftoc3p3 = 100.;
01805
01806 double csVal1 = 0.;
01807 double csVal2 = 0.;
01808
01809 DEBUG_ENTRY( "pah1_fun()" );
01810
01811
01812
01816 x = 1./wavl;
01817
01818 if( x >= x2a )
01819 {
01820
01821 double anu_ev = x/x2a*EVRYD;
01822
01823
01824 t_ADfA::Inst().set_version( PHFIT95 );
01825
01826 term1 = t_ADfA::Inst().phfit(ipHYDROGEN+1,1,1,anu_ev);
01827 term2 = 0.;
01828 for( j=1; j <= 3; ++j )
01829 term2 += t_ADfA::Inst().phfit(ipCARBON+1,6,j,anu_ev);
01830
01831 csVal1 = (xnh*term1 + xnc*term2)*1.e-18;
01832 }
01833
01834 if( x <= 2.*x2a )
01835 {
01836 cval1 = log(sqrt(xnc)*ftoc3p3/1.2328)/12.2;
01837
01838 term = POW2(MIN2(x,x1))*(3.*x1 - 2.*MIN2(x,x1))/POW3(x1);
01839
01840 term1 = (0.1*x + 0.41)*POW2(MAX2(x-x2,0.));
01841
01842
01843
01844
01845
01846
01847
01848
01849
01850
01851
01852
01853
01854
01855
01856
01857
01858
01859
01860
01861
01862
01863
01864 term2 = exp(cval1*(1. - (x1/MIN2(x,x1))));
01865
01866 term3 = p3*exp(-POW2(x/x3))*MIN2(x,x3)/x3;
01867
01868 csVal2 = xnc*((p1*term + p2*term1)*term2 + term3);
01869 }
01870
01871 if( x2a <= x && x <= 2.*x2a )
01872 {
01873
01874 double frac = POW2(2.-x/x2a);
01875 pah1_fun_v = exp((1.-frac)*log(csVal1) + frac*log(csVal2));
01876 }
01877 else
01878 {
01879
01880 pah1_fun_v = csVal1 + csVal2;
01881 }
01882
01883
01884
01885 if( wl1a <= wavl && wl1d >= wavl )
01886 {
01887 if( wavl < wl1b )
01888 term = p4*(wavl - wl1a)/(wl1b - wl1a);
01889 else
01890 term = ( wavl > wl1c ) ? p4*(wl1d - wavl)/(wl1d - wl1c) : p4;
01891 pah1_fun_v += term*xnc;
01892 }
01893 if( wl2a <= wavl && wl2c >= wavl )
01894 {
01895 term = ( wavl < wl2b ) ? p5*((wavl - wl2a)/delwl2) : p5*sqrt(1.-POW2((wavl-wl2a)/(wl2c-wl2a)));
01896 pah1_fun_v += term*xnc;
01897 }
01898 if( wl3m-10.*wl3sig <= wavl && wavl <= wl3m+10.*wl3sig )
01899 {
01900
01901
01902 term = 1.1*pah1_strength[0]*exp(-0.5*POW2((wavl-wl3m)/wl3sig));
01903 pah1_fun_v += term*xnh;
01904 }
01905
01906
01907 for( j=0; j < 7; j++ )
01908 {
01909 term1 = (wavl - pah1_wlBand[j])/pah1_width[j];
01910 term = 0.;
01911 if( j == 2 )
01912 {
01913
01914
01915
01916
01917 if( term1 >= -1. && term1 < -0.5 )
01918 {
01919 term = pah1_strength[j]/(3.*pah1_width[j]);
01920 term *= 2. + 2.*term1;
01921 }
01922 if( term1 >= -0.5 && term1 <= 1.5 )
01923 term = pah1_strength[j]/(3.*pah1_width[j]);
01924 if( term1 > 1.5 && term1 <= 3. )
01925 {
01926 term = pah1_strength[j]/(3.*pah1_width[j]);
01927 term *= 2. - term1*2./3.;
01928 }
01929 }
01930 else
01931 {
01932
01933
01934
01935
01936 if( term1 >= -2. && term1 < -1. )
01937 {
01938 term = pah1_strength[j]/(3.*pah1_width[j]);
01939 term *= 2. + term1;
01940 }
01941 if( term1 >= -1. && term1 <= 1. )
01942 term = pah1_strength[j]/(3.*pah1_width[j]);
01943 if( term1 > 1. && term1 <= 2. )
01944 {
01945 term = pah1_strength[j]/(3.*pah1_width[j]);
01946 term *= 2. - term1;
01947 }
01948 }
01949 if( j == 0 || j > 2 )
01950 term *= xnhoc;
01951 pah1_fun_v += term*xnc;
01952 }
01953
01954 *cs_abs = pah1_fun_v;
01955
01956
01957 *cs_sct = 0.1*pah1_fun_v;
01958 *cosb = 0.;
01959 *error = 0;
01960
01961 DEBUG_EXIT( "pah1_fun()" );
01962
01963 return;
01964 }
01965
01966
01967 static void tbl_fun(double wavl,
01968 sd_data *sd,
01969 grain_data *gd,
01970 double *cs_abs,
01971 double *cs_sct,
01972 double *cosb,
01973 int *error)
01974 {
01975 bool lgOutOfBounds;
01976 long int ind;
01977 double anu = WAVNRYD/wavl*1.e4;
01978
01979 DEBUG_ENTRY( "tbl_fun()" );
01980
01981
01982 if( sd == NULL )
01983 TotalInsanity();
01984
01988 find_arr(anu,gd->opcAnu,gd->nOpcData,&ind,&lgOutOfBounds);
01989 if( !lgOutOfBounds )
01990 {
01991 double a1g;
01992 double frac = log(anu/gd->opcAnu[ind])/log(gd->opcAnu[ind+1]/gd->opcAnu[ind]);
01993
01994 *cs_abs = exp((1.-frac)*log(gd->opcData[0][ind])+frac*log(gd->opcData[0][ind+1]));
01995 ASSERT( *cs_abs > 0. );
01996 if( gd->nOpcCols > 1 )
01997 *cs_sct = exp((1.-frac)*log(gd->opcData[1][ind])+frac*log(gd->opcData[1][ind+1]));
01998 else
01999 *cs_sct = 0.1*(*cs_abs);
02000 ASSERT( *cs_sct > 0. );
02001 if( gd->nOpcCols > 2 )
02002 a1g = exp((1.-frac)*log(gd->opcData[2][ind])+frac*log(gd->opcData[2][ind+1]));
02003 else
02004 a1g = 1.;
02005 ASSERT( a1g > 0. );
02006 *cosb = 1. - a1g;
02007 *error = 0;
02008 }
02009 else
02010 {
02011 *cs_abs = -1.;
02012 *cs_sct = -1.;
02013 *cosb = -2.;
02014 *error = 3;
02015 }
02016
02017 DEBUG_EXIT( "tbl_fun()" );
02018 return;
02019 }
02020
02021
02022 static double size_distr(double size,
02023 sd_data *sd)
02024 {
02025 bool lgOutOfBounds;
02026 long ind;
02027 double frac,
02028 res,
02029 x;
02030
02031 DEBUG_ENTRY( "size_distr()" );
02032
02033 if( size >= sd->lim[ipBLo] && size <= sd->lim[ipBHi] )
02034 switch( sd->sdCase )
02035 {
02036 case SD_SINGLE_SIZE:
02037 res = 1.;
02038 break;
02039 case SD_POWERLAW:
02040
02041 case SD_EXP_CUTOFF1:
02042 case SD_EXP_CUTOFF2:
02043 case SD_EXP_CUTOFF3:
02044
02045
02046 res = pow(size,sd->a[ipExp]);
02047 if( sd->a[ipBeta] < 0. )
02048 res /= (1. - sd->a[ipBeta]*size);
02049 else if( sd->a[ipBeta] > 0. )
02050 res *= (1. + sd->a[ipBeta]*size);
02051 if( size < sd->a[ipBLo] && sd->a[ipSLo] > 0. )
02052 res *= exp(-powi((sd->a[ipBLo]-size)/sd->a[ipSLo],nint(sd->a[ipAlpha])));
02053 if( size > sd->a[ipBHi] && sd->a[ipSHi] > 0. )
02054 res *= exp(-powi((size-sd->a[ipBHi])/sd->a[ipSHi],nint(sd->a[ipAlpha])));
02055 break;
02056 case SD_LOG_NORMAL:
02057 x = log(size/sd->a[ipGCen])/sd->a[ipGSig];
02058 res = exp(-0.5*POW2(x))/size;
02059 break;
02060 case SD_LIN_NORMAL:
02061 x = (size-sd->a[ipGCen])/sd->a[ipGSig];
02062 res = exp(-0.5*POW2(x))/size;
02063 break;
02064 case SD_TABLE:
02065 find_arr(log(size),sd->ln_a,sd->npts,&ind,&lgOutOfBounds);
02066 if( lgOutOfBounds )
02067 {
02068 fprintf( ioQQQ, " size distribution table has insufficient range\n" );
02069 fprintf( ioQQQ, " requested size: %.5f table range %.5f - %.5f\n",
02070 size, exp(sd->ln_a[0]), exp(sd->ln_a[sd->npts-1]) );
02071 puts( "[Stop in size_distr]" );
02072 cdEXIT(EXIT_FAILURE);
02073 }
02074 frac = (log(size)-sd->ln_a[ind])/(sd->ln_a[ind+1]-sd->ln_a[ind]);
02075 ASSERT( frac > 0.-10.*DBL_EPSILON && frac < 1.+10.*DBL_EPSILON );
02076 res = (1.-frac)*sd->ln_a4dNda[ind] + frac*sd->ln_a4dNda[ind+1];
02077
02078 res = exp(res)/POW4(size);
02079 break;
02080 case SD_ILLEGAL:
02081 default:
02082 fprintf( ioQQQ, " insane case for grain size distribution: %d\n" , sd->sdCase );
02083 ShowMe();
02084 puts( "[Stop in size_distr]" );
02085 cdEXIT(EXIT_FAILURE);
02086 }
02087 else
02088 res = 0.;
02089
02090 DEBUG_EXIT( "size_distr()" );
02091 return res;
02092 }
02093
02094
02095
02096
02097
02098
02099 static double search_limit(double ref,
02100 double step,
02101 double rel_cutoff,
02102 sd_data sd)
02103 {
02104 long i;
02105 double f1,
02106 f2,
02107 fmid,
02108 renorm,
02109 x1,
02110 x2 = DBL_MAX,
02111 xmid = DBL_MAX;
02112
02113
02114 const double TOLER = 1.e-6;
02115
02116 DEBUG_ENTRY( "search_limit()" );
02117
02118
02119 ASSERT( rel_cutoff > 0. && rel_cutoff < 1. );
02120
02121 if( step == 0. )
02122 {
02123 DEBUG_EXIT( "search_limit()" );
02124 return ref;
02125 }
02126
02127
02128
02129
02130 sd.lim[ipBLo] = 0.;
02131 sd.lim[ipBHi] = DBL_MAX;
02132
02133 x1 = ref;
02134
02135 f1 = -log(rel_cutoff);
02136 renorm = f1 - log(POW4(x1)*size_distr(x1,&sd));
02137
02138
02139 f2 = 1.;
02140 for( i=0; i < 20 && f2 > 0.; ++i )
02141 {
02142 x2 = MAX2(ref + step,SMALLEST_GRAIN);
02143 f2 = log(POW4(x2)*size_distr(x2,&sd)) + renorm;
02144 if( f2 >= 0. )
02145 {
02146 x1 = x2;
02147 f1 = f2;
02148 }
02149 step *= 2.;
02150 }
02151 if( f2 > 0. )
02152 {
02153 fprintf( ioQQQ, " Could not bracket solution\n" );
02154 puts( "[Stop in search_limit]" );
02155 cdEXIT(EXIT_FAILURE);
02156 }
02157
02158
02159 while( 2.*fabs(x1-x2)/(x1+x2) > TOLER )
02160 {
02161 xmid = (x1+x2)/2.;
02162 fmid = log(POW4(xmid)*size_distr(xmid,&sd)) + renorm;
02163
02164 if( fmid == 0. )
02165 break;
02166
02167 if( f1*fmid > 0. )
02168 {
02169 x1 = xmid;
02170 f1 = fmid;
02171 }
02172 else
02173 {
02174 x2 = xmid;
02175 f2 = fmid;
02176 }
02177 }
02178
02179 DEBUG_EXIT( "search_limit()" );
02180 return (x1+x2)/2.;
02181 }
02182
02183
02184 static void mie_calc_ial( grain_data *gd,
02185 long int n,
02186 double invlen[],
02187 const char *chString,
02188 bool *lgWarning)
02189 {
02190
02191 int *ErrorIndex;
02192 bool lgErrorOccurred=true,
02193 lgOutOfBounds;
02194 long int i,
02195 ind,
02196 j;
02197 double frac,
02198 InvDep,
02199 nim,
02200 wavlen;
02201
02202 DEBUG_ENTRY( "mie_calc_ial()" );
02203
02204
02205 ASSERT( gd->rfiType == RFI_TABLE );
02206
02207
02208 ErrorIndex = (int*)MALLOC(sizeof(int)*(unsigned)rfield.nupper);
02209
02210 for( i=0; i < n; i++ )
02211 {
02212 wavlen = WAVNRYD/rfield.anu[i]*1.e4;
02213
02214 ErrorIndex[i] = 0;
02215 lgErrorOccurred = false;
02216 invlen[i] = 0.;
02217
02218 for( j=0; j < gd->nAxes; j++ )
02219 {
02220
02221 find_arr(wavlen,gd->wavlen[j],gd->ndata[j],&ind,&lgOutOfBounds);
02222 if( lgOutOfBounds )
02223 {
02224 ErrorIndex[i] = 3;
02225 lgErrorOccurred = true;
02226 invlen[i] = 0.;
02227 break;
02228 }
02229 frac = (wavlen-gd->wavlen[j][ind])/(gd->wavlen[j][ind+1]-gd->wavlen[j][ind]);
02230 nim = (1.-frac)*gd->n[j][ind].imag() + frac*gd->n[j][ind+1].imag();
02231
02232
02233 InvDep = PI4*nim/wavlen*1.e4;
02234 ASSERT( InvDep > 0. );
02235
02236 invlen[i] += InvDep*gd->wt[j];
02237 }
02238 }
02239
02240 if( lgErrorOccurred )
02241 {
02242 mie_repair(chString,n,3,3,rfield.anu,invlen,ErrorIndex,false,lgWarning);
02243 }
02244
02245 free( ErrorIndex );
02246
02247 DEBUG_EXIT( "mie_calc_ial()" );
02248 return;
02249 }
02250
02251
02252
02253 #define NPTS_DERIV 8
02254 #define NPTS_COMB (NPTS_DERIV*(NPTS_DERIV-1)/2)
02255
02256
02257 static void mie_repair( const char *chString,
02258 long int n,
02259 int val,
02260 int del,
02261 float anu[],
02262 double data[],
02263 int ErrorIndex[],
02264 bool lgRound,
02265 bool *lgWarning)
02266 {
02267 bool lgExtrapolate,
02268 lgVerbose;
02269 long int i1,
02270 i2,
02271 ind1,
02272 ind2,
02273 j;
02274 double dx,
02275 sgn,
02276 slp1,
02277 xlg1,
02278 xlg2,
02279 y1lg1,
02280 y1lg2;
02281
02282
02283 const long BIG_INTERPOLATION = 10;
02284
02285 DEBUG_ENTRY( "mie_repair()" );
02286
02287 lgVerbose = ( chString[0] != '\0' );
02288
02289 for( ind1=0; ind1 < n; )
02290 {
02291 if( ErrorIndex[ind1] == val )
02292 {
02293
02294 ind2 = ind1;
02295 while( ind2 < n && ErrorIndex[ind2] == val )
02296 ind2++;
02297
02298 if( lgVerbose )
02299 fprintf( ioQQQ, " %s", chString );
02300
02301 if( ind1 == 0 )
02302 {
02303
02304 i1 = ind2;
02305 i2 = ind2+NPTS_DERIV-1;
02306 lgExtrapolate = true;
02307 sgn = +1.;
02308 if( lgVerbose )
02309 {
02310 fprintf( ioQQQ, " extrapolated below %.4e Ryd\n",anu[i1] );
02311 }
02312 }
02313 else if( ind2 == n )
02314 {
02315
02316 i1 = ind1-NPTS_DERIV;
02317 i2 = ind1-1;
02318 lgExtrapolate = true;
02319 sgn = -1.;
02320 if( lgVerbose )
02321 {
02322 fprintf( ioQQQ, " extrapolated above %.4e Ryd\n",anu[i2] );
02323 }
02324 }
02325 else
02326 {
02327
02328 i1 = ind1-1;
02329 i2 = ind2;
02330 lgExtrapolate = false;
02331 sgn = 0.;
02332 if( lgVerbose )
02333 {
02334 fprintf( ioQQQ, " interpolated between %.4e and %.4e Ryd\n",
02335 anu[i1],anu[i2] );
02336 }
02337 if( i2-i1-1 > BIG_INTERPOLATION )
02338 {
02339 if( lgVerbose )
02340 {
02341 fprintf( ioQQQ, " ***Warning: extensive interpolation used\n" );
02342 }
02343 *lgWarning = true;
02344 }
02345 }
02346
02347 if( i1 < 0 || i2 >= n )
02348 {
02349 fprintf( ioQQQ, " Insufficient data for extrapolation\n" );
02350 puts( "[Stop in mie_repair]" );
02351 cdEXIT(EXIT_FAILURE);
02352 }
02353
02354 xlg1 = log(anu[i1]);
02355 y1lg1 = log(data[i1]);
02356
02357 if( lgExtrapolate )
02358 slp1 = mie_find_slope(anu,data,ErrorIndex,i1,i2,val,lgVerbose,lgWarning);
02359 else
02360 {
02361 xlg2 = log(anu[i2]);
02362 y1lg2 = log(data[i2]);
02363 slp1 = (y1lg2-y1lg1)/(xlg2-xlg1);
02364 }
02365 if( lgRound && lgExtrapolate && sgn > 0. )
02366 {
02367
02368
02369
02370 slp1 = MAX2(slp1,0.);
02371 }
02372
02373 else if( lgExtrapolate && sgn*slp1 < 0. )
02374 {
02375 fprintf( ioQQQ, " Illegal value for slope in extrapolation %.6e\n", slp1 );
02376 puts( "[Stop in mie_repair]" );
02377 cdEXIT(EXIT_FAILURE);
02378 }
02379
02380 for( j=ind1; j < ind2; j++ )
02381 {
02382 dx = log(anu[j]) - xlg1;
02383 data[j] = exp(y1lg1 + dx*slp1);
02384 ErrorIndex[j] -= del;
02385 }
02386
02387 ind1 = ind2;
02388 }
02389 else
02390 {
02391 ind1++;
02392 }
02393 }
02394
02395 for( j=0; j < n; j++ )
02396 {
02397 if( ErrorIndex[j] > val-del )
02398 {
02399 fprintf( ioQQQ, " Internal error in mie_repair, index=%ld, val=%d\n",j,ErrorIndex[j] );
02400 ShowMe();
02401 puts( "[Stop in mie_repair]" );
02402 cdEXIT(EXIT_FAILURE);
02403 }
02404 }
02405
02406 DEBUG_EXIT( "mie_repair()" );
02407 return;
02408 }
02409
02410
02411 static double mie_find_slope( const float anu[],
02412 const double data[],
02413 const int ErrorIndex[],
02414 long i1,
02415 long i2,
02416 int val,
02417 bool lgVerbose,
02418 bool *lgWarning)
02419 {
02420 long i,
02421 j,
02422 k;
02423 double s1,
02424 s2,
02425 slope,
02426 slp1[NPTS_COMB],
02427 stdev;
02428
02429
02430
02431 const double LARGE_STDEV = 0.2;
02432
02433 DEBUG_ENTRY( "mie_find_slope()" );
02434
02435
02436 ASSERT( i2-i1 == NPTS_DERIV-1 );
02437 for( i=i1; i <= i2; i++ )
02438 {
02439 ASSERT( ErrorIndex[i] < val );
02440 ASSERT( anu[i] > 0. && data[i] > 0. );
02441 }
02442
02443
02444 for( i=0; i < NPTS_COMB; i++ )
02445 {
02446 slp1[i] = -DBL_MAX;
02447 }
02448
02449 k = 0;
02450
02451
02452 for( i=i1; i < i2; i++ )
02453 {
02454 for( j=i+1; j <= i2; j++ )
02455 {
02456 slp1[k++] = log(data[j]/data[i])/log(anu[j]/anu[i]);
02457 }
02458 }
02459
02460 for( i=0; i <= NPTS_COMB/2; i++ )
02461 {
02462 for( j=i+1; j < NPTS_COMB; j++ )
02463 {
02464 if( slp1[i] > slp1[j] )
02465 {
02466 double xxx = slp1[i];
02467 slp1[i] = slp1[j];
02468 slp1[j] = xxx;
02469 }
02470 }
02471 }
02472
02473 slope = ( NPTS_COMB%2 == 1 ) ? slp1[NPTS_COMB/2] : (slp1[NPTS_COMB/2-1] + slp1[NPTS_COMB/2])/2.;
02474
02475
02476 s1 = s2 = 0.;
02477 for( i=0; i < NPTS_COMB; i++ )
02478 {
02479 s1 += slp1[i];
02480 s2 += POW2(slp1[i]);
02481 }
02482
02483 stdev = MAX2(s2/(double)NPTS_COMB - POW2(s1/(double)NPTS_COMB),0.);
02484 stdev = sqrt(stdev);
02485
02486 #if 0
02487 for( i=i1; i <= i2; i++ )
02488 printf("input: %ld %.4e %.4e\n",i,anu[i],data[i]);
02489 for( i=0; i < NPTS_COMB; i++ )
02490 printf("%.3f ",slp1[i]);
02491 printf("\n");
02492 printf("slope %.3f +/- %.3f\n",slope,stdev);
02493 #endif
02494
02495
02496 if( stdev > LARGE_STDEV )
02497 {
02498 if( lgVerbose )
02499 {
02500 fprintf( ioQQQ, " ***Warning: slope for extrapolation may be unreliable\n" );
02501 }
02502 *lgWarning = true;
02503 }
02504
02505 DEBUG_EXIT( "mie_find_slope()" );
02506 return slope;
02507 }
02508
02509
02510
02511 static void mie_read_rfi( const char *chFile,
02512 grain_data *gd)
02513 {
02514 bool lgLogData=false;
02515 long int dl,
02516 help,
02517 i,
02518 nelem,
02519 j,
02520 nridf,
02521 sgn = 0;
02522 double eps1,
02523 eps2,
02524 LargestLog,
02525 molw,
02526 nAtoms,
02527 nr,
02528 ni,
02529 tmp1,
02530 tmp2,
02531 total = 0.;
02532 char chLine[FILENAME_PATH_LENGTH_2],
02533 chWord[FILENAME_PATH_LENGTH_2];
02534 FILE *io2;
02535
02536 DEBUG_ENTRY( "mie_read_rfi()" );
02537
02538 if( (io2 = fopen(chFile,"r")) == NULL )
02539 {
02540 fprintf( ioQQQ, " Could not open %s for reading\n",chFile);
02541 puts( "[Stop in mie_read_rfi]" );
02542 cdEXIT(EXIT_FAILURE);
02543 }
02544
02545 dl = 0;
02546
02547
02548 mie_next_data(chFile,io2,chLine,&dl);
02549 mie_read_long(chFile,chLine,&gd->magic,true,dl);
02550 if( gd->magic != MAGIC_RFI )
02551 {
02552 fprintf( ioQQQ, " Refractive index file %s has obsolete magic number\n",chFile );
02553 fprintf( ioQQQ, " I found %ld, but expected %ld on line #%ld\n",gd->magic,MAGIC_RFI,dl );
02554 fprintf( ioQQQ, " Please replace this file with an up to date version\n" );
02555 puts( "[Stop in mie_read_rfi]" );
02556 cdEXIT(EXIT_FAILURE);
02557 }
02558
02559
02560 mie_next_data(chFile,io2,chLine,&dl);
02561 mie_read_word(chLine,chWord,FILENAME_PATH_LENGTH_2,false);
02562 mie_read_form(chWord,gd->elmAbun,&nAtoms,&molw);
02563
02564
02565 gd->mol_weight = molw;
02566 gd->atom_weight = gd->mol_weight/nAtoms;
02567
02568
02569 mie_next_data(chFile,io2,chLine,&dl);
02570 mie_read_double(chFile,chLine,&gd->abun,true,dl);
02571
02572
02573 mie_next_data(chFile,io2,chLine,&dl);
02574 mie_read_double(chFile,chLine,&gd->depl,true,dl);
02575 if( gd->depl > 1. )
02576 {
02577 fprintf( ioQQQ, " Illegal value for default depletion in %s\n",chFile );
02578 fprintf( ioQQQ, " Line #%ld, depl=%14.6e\n",dl,gd->depl);
02579 puts( "[Stop in mie_read_rfi]" );
02580 cdEXIT(EXIT_FAILURE);
02581 }
02582
02583 for( nelem=0; nelem < LIMELM; nelem++ )
02584 gd->elmAbun[nelem] *= gd->abun*gd->depl;
02585
02586
02587 mie_next_data(chFile,io2,chLine,&dl);
02588 mie_read_double(chFile,chLine,&gd->rho,false,dl);
02589
02590
02591 mie_next_data(chFile,io2,chLine,&dl);
02592 mie_read_long(chFile,chLine,&help,true,dl);
02593 gd->matType = (mat_type)help;
02594 if( gd->matType >= MAT_TOP )
02595 {
02596 fprintf( ioQQQ, " Illegal value for material type in %s\n",chFile );
02597 fprintf( ioQQQ, " Line #%ld, type=%d\n",dl,gd->matType);
02598 puts( "[Stop in mie_read_rfi]" );
02599 cdEXIT(EXIT_FAILURE);
02600 }
02601
02602
02603 mie_next_data(chFile,io2,chLine,&dl);
02604 mie_read_double(chFile,chLine,&gd->work,true,dl);
02605
02606
02607 mie_next_data(chFile,io2,chLine,&dl);
02608 mie_read_double(chFile,chLine,&gd->bandgap,false,dl);
02609 if( gd->bandgap >= gd->work )
02610 {
02611 fprintf( ioQQQ, " Illegal value for bandgap in %s\n",chFile );
02612 fprintf( ioQQQ, " Line #%ld, bandgap=%.4e, work function=%.4e\n",dl,gd->bandgap,gd->work);
02613 fprintf( ioQQQ, " Bandgap should always be less than work function\n");
02614 puts( "[Stop in mie_read_rfi]" );
02615 cdEXIT(EXIT_FAILURE);
02616 }
02617
02618
02619 mie_next_data(chFile,io2,chLine,&dl);
02620 mie_read_double(chFile,chLine,&gd->therm_eff,true,dl);
02621 if( gd->therm_eff > 1.f )
02622 {
02623 fprintf( ioQQQ, " Illegal value for thermionic efficiency in %s\n",chFile );
02624 fprintf( ioQQQ, " Line #%ld, value=%.4e\n",dl,gd->therm_eff);
02625 fprintf( ioQQQ, " Allowed values are 0. < efficiency <= 1.\n");
02626 puts( "[Stop in mie_read_rfi]" );
02627 cdEXIT(EXIT_FAILURE);
02628 }
02629
02630
02631 mie_next_data(chFile,io2,chLine,&dl);
02632 mie_read_double(chFile,chLine,&gd->subl_temp,true,dl);
02633
02634
02635 mie_next_data(chFile,io2,chLine,&dl);
02636 mie_read_word(chLine,chWord,WORDLEN,true);
02637
02638 if( nMatch( "RFI_", chWord ) )
02639 gd->rfiType = RFI_TABLE;
02640 else if( nMatch( "OPC_", chWord ) )
02641 gd->rfiType = OPC_TABLE;
02642 else if( nMatch( "GREY", chWord ) )
02643 gd->rfiType = OPC_GREY;
02644 else if( nMatch( "PAH1", chWord ) )
02645 gd->rfiType = OPC_PAH1;
02646 else
02647 {
02648 fprintf( ioQQQ, " Illegal keyword in %s\n",chFile );
02649 fprintf( ioQQQ, " Line #%ld, value=%s\n",dl,chWord);
02650 fprintf( ioQQQ, " Allowed values are: RFI_TBL, OPC_TBL, GREY, PAH1\n");
02651 puts( "[Stop in mie_read_rfi]" );
02652 cdEXIT(EXIT_FAILURE);
02653 }
02654
02655 switch( gd->rfiType )
02656 {
02657 case RFI_TABLE:
02658
02659
02660 mie_next_data(chFile,io2,chLine,&dl);
02661 mie_read_long(chFile,chLine,&nridf,true,dl);
02662 if( nridf > 3 )
02663 {
02664 fprintf( ioQQQ, " Illegal data code in %s\n",chFile );
02665 fprintf( ioQQQ, " Line #%ld, data code=%ld\n",dl,nridf);
02666 puts( "[Stop in mie_read_rfi]" );
02667 cdEXIT(EXIT_FAILURE);
02668 }
02669
02670
02671
02672 mie_next_data(chFile,io2,chLine,&dl);
02673 mie_read_long(chFile,chLine,&gd->nAxes,true,dl);
02674 if( gd->nAxes > NAX )
02675 {
02676 fprintf( ioQQQ, " Illegal no. of axes in %s\n",chFile );
02677 fprintf( ioQQQ, " Line #%ld, number=%ld\n",dl,gd->nAxes);
02678 puts( "[Stop in mie_read_rfi]" );
02679 cdEXIT(EXIT_FAILURE);
02680 }
02681
02682
02683 mie_next_data(chFile,io2,chLine,&dl);
02684 switch( gd->nAxes )
02685 {
02686 case 1:
02687 mie_read_double(chFile,chLine,&gd->wt[0],true,dl);
02688 total = gd->wt[0];
02689 break;
02690 case 2:
02691 if( sscanf( chLine, "%lf %lf", &gd->wt[0], &gd->wt[1] ) != 2 )
02692 {
02693 fprintf( ioQQQ, " Syntax error in %s\n",chFile);
02694 fprintf( ioQQQ, " Line #%ld: %s\n",dl,chLine);
02695 puts( "[Stop in mie_read_rfi]" );
02696 cdEXIT(EXIT_FAILURE);
02697 }
02698 if( gd->wt[0] <= 0. || gd->wt[1] <= 0. )
02699 {
02700 fprintf( ioQQQ, " Illegal data in %s\n",chFile);
02701 fprintf( ioQQQ, " Line #%ld: %s\n",dl,chLine);
02702 puts( "[Stop in mie_read_rfi]" );
02703 cdEXIT(EXIT_FAILURE);
02704 }
02705 total = gd->wt[0] + gd->wt[1];
02706 break;
02707 case 3:
02708 if( sscanf( chLine, "%lf %lf %lf", &gd->wt[0], &gd->wt[1], &gd->wt[2] ) != 3 )
02709 {
02710 fprintf( ioQQQ, " Syntax error in %s\n",chFile);
02711 fprintf( ioQQQ, " Line #%ld: %s\n",dl,chLine);
02712 puts( "[Stop in mie_read_rfi]" );
02713 cdEXIT(EXIT_FAILURE);
02714 }
02715 if( gd->wt[0] <= 0. || gd->wt[1] <= 0. || gd->wt[2] <= 0. )
02716 {
02717 fprintf( ioQQQ, " Illegal data in %s\n",chFile);
02718 fprintf( ioQQQ, " Line #%ld: %s\n",dl,chLine);
02719 puts( "[Stop in mie_read_rfi]" );
02720 cdEXIT(EXIT_FAILURE);
02721 }
02722 total = gd->wt[0] + gd->wt[1] + gd->wt[2];
02723 break;
02724 default:
02725 fprintf( ioQQQ, " insane case for gd->nAxes: %ld\n", gd->nAxes );
02726 ShowMe();
02727 puts( "[Stop in mie_read_rfi]" );
02728 cdEXIT(EXIT_FAILURE);
02729 }
02730 for( j=0; j < gd->nAxes; j++ )
02731 {
02732 gd->wt[j] /= total;
02733
02734
02735 mie_next_data(chFile,io2,chLine,&dl);
02736 mie_read_long(chFile,chLine,&gd->ndata[j],false,dl);
02737 if( gd->ndata[j] < 2 )
02738 {
02739 fprintf( ioQQQ, " Illegal number of data points in %s\n",chFile );
02740 fprintf( ioQQQ, " Line #%ld, number=%ld\n",dl,gd->ndata[j]);
02741 puts( "[Stop in mie_read_rfi]" );
02742 cdEXIT(EXIT_FAILURE);
02743 }
02744
02745
02746
02747 gd->wavlen[j] = (double*)MALLOC(sizeof(double)*(unsigned)gd->ndata[j]);
02748 gd->n[j] = (complex<double>*)MALLOC(sizeof(complex<double>)*(unsigned)gd->ndata[j]);
02749 gd->nr1[j] = (double*)MALLOC(sizeof(double)*(unsigned)gd->ndata[j]);
02750
02751 for( i=0; i < gd->ndata[j]; i++ )
02752 {
02753
02754
02755 mie_next_data(chFile,io2,chLine,&dl);
02756 if( sscanf( chLine, "%lf %lf %lf", gd->wavlen[j]+i, &nr, &ni ) != 3 )
02757 {
02758 fprintf( ioQQQ, " Syntax error in %s\n",chFile);
02759 fprintf( ioQQQ, " Line #%ld: %s\n",dl,chLine);
02760 puts( "[Stop in mie_read_rfi]" );
02761 cdEXIT(EXIT_FAILURE);
02762 }
02763 if( gd->wavlen[j][i] <= 0. )
02764 {
02765 fprintf( ioQQQ, " Illegal value for wavelength in %s\n",chFile);
02766 fprintf( ioQQQ, " Line #%ld, wavl=%14.6e\n",dl,gd->wavlen[j][i]);
02767 puts( "[Stop in mie_read_rfi]" );
02768 cdEXIT(EXIT_FAILURE);
02769 }
02770
02771
02772 if( i == 1 )
02773 {
02774 sgn = sign3(gd->wavlen[j][1]-gd->wavlen[j][0]);
02775 if( sgn == 0 )
02776 {
02777 fprintf( ioQQQ, " Illegal value for wavelength in %s\n",chFile);
02778 fprintf( ioQQQ, " Line #%ld, wavl=%14.6e\n",dl,gd->wavlen[j][i]);
02779 puts( "[Stop in mie_read_rfi]" );
02780 cdEXIT(EXIT_FAILURE);
02781 }
02782 }
02783 else if( i > 1 )
02784 {
02785 if( sign3(gd->wavlen[j][i]-gd->wavlen[j][i-1]) != sgn )
02786 {
02787 fprintf( ioQQQ, " Illegal value for wavelength in %s\n",chFile);
02788 fprintf( ioQQQ, " Line #%ld, wavl=%14.6e\n",dl,gd->wavlen[j][i]);
02789 puts( "[Stop in mie_read_rfi]" );
02790 cdEXIT(EXIT_FAILURE);
02791 }
02792 }
02793
02794
02795
02796 switch( nridf )
02797 {
02798 case 1:
02799 eps1 = nr;
02800 eps2 = ni;
02801 dftori(&nr,&ni,eps1,eps2);
02802 gd->nr1[j][i] = nr - 1.;
02803 break;
02804 case 2:
02805 gd->nr1[j][i] = nr;
02806 nr += 1.;
02807 break;
02808 case 3:
02809 gd->nr1[j][i] = nr - 1.;
02810 break;
02811 default:
02812 fprintf( ioQQQ, " insane case for nridf: %ld\n", nridf );
02813 ShowMe();
02814 puts( "[Stop in mie_read_rfi]" );
02815 cdEXIT(EXIT_FAILURE);
02816 }
02817 gd->n[j][i] = complex<double>(nr,ni);
02818
02819
02820 if( nr <= 0. || ni < 0. )
02821 {
02822 fprintf( ioQQQ, " Illegal value for refractive index in %s\n",chFile);
02823 fprintf( ioQQQ, " Line #%ld, (nr,ni)=(%14.6e,%14.6e)\n",dl,nr,ni);
02824 puts( "[Stop in mie_read_rfi]" );
02825 cdEXIT(EXIT_FAILURE);
02826 }
02827 ASSERT( fabs(nr-1.-gd->nr1[j][i]) < 10.*nr*DBL_EPSILON );
02828 }
02829 }
02830 break;
02831 case OPC_TABLE:
02832
02833
02834
02835
02836
02837
02838 mie_next_data(chFile,io2,chLine,&dl);
02839 mie_read_long(chFile,chLine,&gd->nOpcCols,true,dl);
02840 if( gd->nOpcCols > NDAT )
02841 {
02842 fprintf( ioQQQ, " Illegal no. of data columns in %s\n",chFile );
02843 fprintf( ioQQQ, " Line #%ld, number=%ld\n",dl,gd->nOpcCols);
02844 puts( "[Stop in mie_read_rfi]" );
02845 cdEXIT(EXIT_FAILURE);
02846 }
02847
02848
02849 mie_next_data(chFile,io2,chLine,&dl);
02850 mie_read_word(chLine,chWord,WORDLEN,true);
02851
02852 if( nMatch( "LINE", chWord ) )
02853 lgLogData = false;
02854 else if( nMatch( "LOG", chWord ) )
02855 lgLogData = true;
02856 else
02857 {
02858 fprintf( ioQQQ, " Keyword not recognized in %s\n",chFile );
02859 fprintf( ioQQQ, " Line #%ld, keyword=%s\n",dl,chWord);
02860 puts( "[Stop in mie_read_rfi]" );
02861 cdEXIT(EXIT_FAILURE);
02862 }
02863
02864
02865
02866 mie_next_data(chFile,io2,chLine,&dl);
02867 mie_read_long(chFile,chLine,&gd->nOpcData,false,dl);
02868 if( gd->nOpcData < 2 )
02869 {
02870 fprintf( ioQQQ, " Illegal number of data points in %s\n",chFile );
02871 fprintf( ioQQQ, " Line #%ld, number=%ld\n",dl,gd->nOpcData);
02872 puts( "[Stop in mie_read_rfi]" );
02873 cdEXIT(EXIT_FAILURE);
02874 }
02875
02876
02877
02878 gd->opcAnu = (double*)MALLOC(sizeof(double)*(unsigned)gd->nOpcData);
02879 for( j=0; j < gd->nOpcCols; j++ )
02880 {
02881 gd->opcData[j] = (double*)MALLOC(sizeof(double)*(unsigned)gd->nOpcData);
02882 }
02883
02884 tmp1 = -log10(1.01*DBL_MIN);
02885 tmp2 = log10(0.99*DBL_MAX);
02886 LargestLog = MIN2(tmp1,tmp2);
02887
02888
02889
02890
02891
02892
02893
02894
02895
02896
02897 for( i=0; i < gd->nOpcData; i++ )
02898 {
02899 mie_next_data(chFile,io2,chLine,&dl);
02900 switch( gd->nOpcCols )
02901 {
02902 case 1:
02903 if( sscanf( chLine, "%lf %lf", &gd->opcAnu[i], &gd->opcData[0][i] ) != 2 )
02904 {
02905 fprintf( ioQQQ, " Syntax error in %s\n",chFile);
02906 fprintf( ioQQQ, " Line #%ld: %s\n",dl,chLine);
02907 puts( "[Stop in mie_read_rfi]" );
02908 cdEXIT(EXIT_FAILURE);
02909 }
02910 break;
02911 case 2:
02912 if( sscanf( chLine, "%lf %lf %lf", &gd->opcAnu[i], &gd->opcData[0][i],
02913 &gd->opcData[1][i] ) != 3 )
02914 {
02915 fprintf( ioQQQ, " Syntax error in %s\n",chFile);
02916 fprintf( ioQQQ, " Line #%ld: %s\n",dl,chLine);
02917 puts( "[Stop in mie_read_rfi]" );
02918 cdEXIT(EXIT_FAILURE);
02919 }
02920 break;
02921 case 3:
02922 if( sscanf( chLine, "%lf %lf %lf %lf", &gd->opcAnu[i], &gd->opcData[0][i],
02923 &gd->opcData[1][i], &gd->opcData[2][i] ) != 4 )
02924 {
02925 fprintf( ioQQQ, " Syntax error in %s\n",chFile);
02926 fprintf( ioQQQ, " Line #%ld: %s\n",dl,chLine);
02927 puts( "[Stop in mie_read_rfi]" );
02928 cdEXIT(EXIT_FAILURE);
02929 }
02930 break;
02931 case 4:
02932 if( sscanf( chLine, "%lf %lf %lf %lf %lf", &gd->opcAnu[i], &gd->opcData[0][i],
02933 &gd->opcData[1][i], &gd->opcData[2][i], &gd->opcData[3][i] ) != 5 )
02934 {
02935 fprintf( ioQQQ, " Syntax error in %s\n",chFile);
02936 fprintf( ioQQQ, " Line #%ld: %s\n",dl,chLine);
02937 puts( "[Stop in mie_read_rfi]" );
02938 cdEXIT(EXIT_FAILURE);
02939 }
02940 break;
02941 default:
02942 fprintf( ioQQQ, " insane case for gd->nOpcCols: %ld\n", gd->nOpcCols );
02943 ShowMe();
02944 puts( "[Stop in mie_read_rfi]" );
02945 cdEXIT(EXIT_FAILURE);
02946 }
02947 if( lgLogData )
02948 {
02949
02950 if( fabs(gd->opcAnu[i]) > LargestLog )
02951 {
02952 fprintf( ioQQQ, " Illegal value for frequency in %s\n",chFile );
02953 fprintf( ioQQQ, " Line #%ld, freq=%14.6e\n",dl,gd->opcAnu[i] );
02954 puts( "[Stop in mie_read_rfi]" );
02955 cdEXIT(EXIT_FAILURE);
02956 }
02957 gd->opcAnu[i] = pow(10.,gd->opcAnu[i]);
02958 for( j=0; j < gd->nOpcCols; j++ )
02959 {
02960 if( fabs(gd->opcData[j][i]) > LargestLog )
02961 {
02962 fprintf( ioQQQ, " Illegal data value in %s\n",chFile );
02963 fprintf( ioQQQ, " Line #%ld, value=%14.6e\n",dl,gd->opcData[j][i] );
02964 puts( "[Stop in mie_read_rfi]" );
02965 cdEXIT(EXIT_FAILURE);
02966 }
02967 gd->opcData[j][i] = pow(10.,gd->opcData[j][i]);
02968 }
02969 }
02970 if( gd->opcAnu[i] <= 0. )
02971 {
02972 fprintf( ioQQQ, " Illegal value for frequency in %s\n",chFile );
02973 fprintf( ioQQQ, " Line #%ld, freq=%14.6e\n",dl,gd->opcAnu[i] );
02974 puts( "[Stop in mie_read_rfi]" );
02975 cdEXIT(EXIT_FAILURE);
02976 }
02977 for( j=0; j < gd->nOpcCols; j++ )
02978 {
02979 if( gd->opcData[j][i] <= 0. )
02980 {
02981 fprintf( ioQQQ, " Illegal data value in %s\n",chFile );
02982 fprintf( ioQQQ, " Line #%ld, value=%14.6e\n",dl,gd->opcData[j][i] );
02983 puts( "[Stop in mie_read_rfi]" );
02984 cdEXIT(EXIT_FAILURE);
02985 }
02986 }
02987
02988
02989 if( i == 1 )
02990 {
02991 sgn = sign3(gd->opcAnu[1]-gd->opcAnu[0]);
02992 if( sgn == 0 )
02993 {
02994 double dataVal = lgLogData ? log10(gd->opcAnu[1]) : gd->opcAnu[1];
02995 fprintf( ioQQQ, " Illegal value for frequency in %s\n",chFile );
02996 fprintf( ioQQQ, " Line #%ld, freq=%14.6e\n",dl,dataVal );
02997 puts( "[Stop in mie_read_rfi]" );
02998 cdEXIT(EXIT_FAILURE);
02999 }
03000 }
03001 else if( i > 1 )
03002 {
03003 if( sign3(gd->opcAnu[i]-gd->opcAnu[i-1]) != sgn )
03004 {
03005 double dataVal = lgLogData ? log10(gd->opcAnu[i]) : gd->opcAnu[i];
03006 fprintf( ioQQQ, " Illegal value for frequency in %s\n",chFile);
03007 fprintf( ioQQQ, " Line #%ld, freq=%14.6e\n",dl,dataVal);
03008 puts( "[Stop in mie_read_rfi]" );
03009 cdEXIT(EXIT_FAILURE);
03010 }
03011 }
03012 }
03013 gd->nAxes = 1;
03014 break;
03015 case OPC_GREY:
03016 case OPC_PAH1:
03017
03018
03019 gd->nAxes = 1;
03020 break;
03021 default:
03022 fprintf( ioQQQ, " Insane value for gd->rfiType: %d, bailing out\n", gd->rfiType );
03023 ShowMe();
03024 puts( "[Stop in mie_read_rfi]" );
03025 cdEXIT(EXIT_FAILURE);
03026 }
03027
03028 fclose(io2);
03029
03030 DEBUG_EXIT( "mie_read_rfi()" );
03031 return;
03032 }
03033
03034
03035
03036 static void mie_read_mix( const char *chFile,
03037 grain_data *gd)
03038 {
03039 emt_type EMTtype;
03040 long int dl,
03041 i,
03042 j,
03043 k,
03044 l,
03045 nelem,
03046 nMaterial,
03047 sumAxes;
03048 double *delta,
03049 *frac,
03050 *frac2,
03051 *frdelta,
03052 maxIndex = DBL_MAX,
03053 minIndex = DBL_MAX,
03054 nAtoms,
03055 sum,
03056 sum2,
03057 wavHi,
03058 wavLo,
03059 wavStep;
03060 complex<double> *eps,
03061 eps_eff(-DBL_MAX,-DBL_MAX);
03062 char chLine[FILENAME_PATH_LENGTH_2],
03063 chWord[FILENAME_PATH_LENGTH_2],
03064 *str;
03065 grain_data *gdArr;
03066 FILE *io2;
03067
03068 DEBUG_ENTRY( "mie_read_mix()" );
03069
03070 if( (io2 = fopen(chFile,"r")) == NULL )
03071 {
03072 fprintf( ioQQQ, " Could not open %s for reading\n",chFile);
03073 puts( "[Stop in mie_read_mix]" );
03074 cdEXIT(EXIT_FAILURE);
03075 }
03076
03077 dl = 0;
03078
03079
03080 mie_next_data(chFile,io2,chLine,&dl);
03081 mie_read_long(chFile,chLine,&gd->magic,true,dl);
03082 if( gd->magic != MAGIC_MIX )
03083 {
03084 fprintf( ioQQQ, " Mixed grain file %s has obsolete magic number\n",chFile );
03085 fprintf( ioQQQ, " I found %ld, but expected %ld on line #%ld\n",gd->magic,MAGIC_MIX,dl );
03086 fprintf( ioQQQ, " Please replace this file with an up to date version\n" );
03087 puts( "[Stop in mie_read_mix]" );
03088 cdEXIT(EXIT_FAILURE);
03089 }
03090
03091
03092 mie_next_data(chFile,io2,chLine,&dl);
03093 mie_read_double(chFile,chLine,&gd->depl,true,dl);
03094 if( gd->depl > 1. )
03095 {
03096 fprintf( ioQQQ, " Illegal value for default depletion in %s\n",chFile );
03097 fprintf( ioQQQ, " Line #%ld, depl=%14.6e\n",dl,gd->depl);
03098 puts( "[Stop in mie_read_mix]" );
03099 cdEXIT(EXIT_FAILURE);
03100 }
03101
03102
03103 mie_next_data(chFile,io2,chLine,&dl);
03104 mie_read_long(chFile,chLine,&nMaterial,true,dl);
03105 if( nMaterial < 2 )
03106 {
03107 fprintf( ioQQQ, " Illegal number of materials in mixed grain file %s\n",chFile );
03108 fprintf( ioQQQ, " I found %ld on line #%ld\n",nMaterial,dl );
03109 fprintf( ioQQQ, " This number should be at least 2\n" );
03110 puts( "[Stop in mie_read_mix]" );
03111 cdEXIT(EXIT_FAILURE);
03112 }
03113
03114 frac = (double*)MALLOC(sizeof(double)*(unsigned)nMaterial);
03115 frac2 = (double*)MALLOC(sizeof(double)*(unsigned)nMaterial);
03116 gdArr = (grain_data*)MALLOC(sizeof(grain_data)*(unsigned)nMaterial);
03117
03118 sum = 0.;
03119 sum2 = 0.;
03120 sumAxes = 0;
03121 for( i=0; i < nMaterial; i++ )
03122 {
03123 char chFile2[FILENAME_PATH_LENGTH_2];
03124
03125 gdArr[i].nAxes = 0;
03126 for( j=0; j < NAX; j++ )
03127 {
03128 gdArr[i].wavlen[j] = NULL;
03129 gdArr[i].n[j] = NULL;
03130 gdArr[i].nr1[j] = NULL;
03131 }
03132 gdArr[i].nOpcCols = 0;
03133 gdArr[i].opcAnu = NULL;
03134 for( j=0; j < NDAT; j++ )
03135 {
03136 gdArr[i].opcData[j] = NULL;
03137 }
03138
03139
03140
03141 mie_next_data(chFile,io2,chLine,&dl);
03142 mie_read_double(chFile,chLine,&frac[i],true,dl);
03143
03144 sum += frac[i];
03145
03146 str = strchr( chLine, '\"' );
03147 if( str != NULL )
03148 {
03149 strcpy( chFile2, ++str );
03150 str = strchr( chFile2, '\"' );
03151 if( str != NULL )
03152 *str = '\0';
03153 }
03154 if( str == NULL )
03155 {
03156 fprintf( ioQQQ, " No pair of double quotes was found on line #%ld of file %s\n",dl,chFile );
03157 fprintf( ioQQQ, " Please supply the refractive index file name between double quotes\n" );
03158 puts( "[Stop in mie_read_mix]" );
03159 cdEXIT(EXIT_FAILURE);
03160 }
03161
03162 mie_read_rfi( chFile2, &gdArr[i] );
03163 if( gdArr[i].rfiType != RFI_TABLE )
03164 {
03165 fprintf( ioQQQ, " Input error on line #%ld of file %s\n",dl,chFile );
03166 fprintf( ioQQQ, " File %s is not of type RFI_TBL, this is illegal\n",chFile2 );
03167 puts( "[Stop in mie_read_mix]" );
03168 cdEXIT(EXIT_FAILURE);
03169 }
03170
03171
03172 if( i == nMaterial-1 && gdArr[i].mol_weight == 0. )
03173 {
03174 fprintf( ioQQQ, " Input error on line #%ld of file %s\n",dl,chFile );
03175 fprintf( ioQQQ, " Last entry in list of materials is vacuum, this is not allowed\n" );
03176 fprintf( ioQQQ, " Please move this entry to an earlier position in the file\n" );
03177 puts( "[Stop in mie_read_mix]" );
03178 cdEXIT(EXIT_FAILURE);
03179 }
03180
03181 frac2[i] = ( gdArr[i].mol_weight > 0. ) ? frac[i] : 0.;
03182 sum2 += frac2[i];
03183
03184 sumAxes += gdArr[i].nAxes;
03185 }
03186
03187
03188 mie_next_data(chFile,io2,chLine,&dl);
03189 mie_read_word(chLine,chWord,WORDLEN,true);
03190
03191 if( nMatch( "FA00", chWord ) )
03192 EMTtype = FARAFONOV00;
03193 else if( nMatch( "ST95", chWord ) )
03194 EMTtype = STOGNIENKO95;
03195 else if( nMatch( "BR35", chWord ) )
03196 EMTtype = BRUGGEMAN35;
03197 else
03198 {
03199 fprintf( ioQQQ, " Keyword not recognized in %s\n",chFile );
03200 fprintf( ioQQQ, " Line #%ld, keyword=%s\n",dl,chWord);
03201 puts( "[Stop in mie_read_mix]" );
03202 cdEXIT(EXIT_FAILURE);
03203 }
03204
03205
03206 for( i=0; i < nMaterial; i++ )
03207 {
03208 frac[i] /= sum;
03209 frac2[i] /= sum2;
03210
03211 for( nelem=0; nelem < LIMELM; nelem++ )
03212 {
03213 gdArr[i].elmAbun[nelem] /= gdArr[i].abun*gdArr[i].depl;
03214 }
03215 }
03216
03217 wavLo = 0.;
03218 wavHi = DBL_MAX;
03219 gd->abun = DBL_MAX;
03220 for( nelem=0; nelem < LIMELM; nelem++ )
03221 {
03222 gd->elmAbun[nelem] = 0.;
03223 }
03224 gd->mol_weight = 0.;
03225 gd->rho = 0.;
03226 gd->work = DBL_MAX;
03227 gd->bandgap = DBL_MAX;
03228 gd->therm_eff = 0.;
03229 gd->subl_temp = DBL_MAX;
03230 gd->nAxes = 1;
03231 gd->wt[0] = 1.;
03232 gd->ndata[0] = MIX_TABLE_SIZE;
03233 gd->rfiType = RFI_TABLE;
03234
03235 for( i=0; i < nMaterial; i++ )
03236 {
03237 for( k=0; k < gdArr[i].nAxes; k++ )
03238 {
03239 double wavMin = MIN2(gdArr[i].wavlen[k][0],gdArr[i].wavlen[k][gdArr[i].ndata[k]-1]);
03240 double wavMax = MAX2(gdArr[i].wavlen[k][0],gdArr[i].wavlen[k][gdArr[i].ndata[k]-1]);
03241 wavLo = MAX2(wavLo,wavMin);
03242 wavHi = MIN2(wavHi,wavMax);
03243 }
03244 minIndex = DBL_MAX;
03245 maxIndex = 0.;
03246 for( nelem=0; nelem < LIMELM; nelem++ )
03247 {
03248 gd->elmAbun[nelem] += frac[i]*gdArr[i].elmAbun[nelem];
03249 if( gd->elmAbun[nelem] > 0. )
03250 {
03251 minIndex = MIN2(minIndex,gd->elmAbun[nelem]);
03252 }
03253 maxIndex = MAX2(maxIndex,gd->elmAbun[nelem]);
03254 }
03255 gd->mol_weight += frac[i]*gdArr[i].mol_weight;
03256 gd->rho += frac[i]*gdArr[i].rho;
03257
03258 if( gdArr[i].mol_weight > 0. )
03259 {
03260 gd->abun = MIN2(gd->abun,gdArr[i].abun/frac[i]);
03261 switch( EMTtype )
03262 {
03263 case FARAFONOV00:
03264
03265 gd->work = gdArr[i].work;
03266 gd->bandgap = gdArr[i].bandgap;
03267 gd->therm_eff = gdArr[i].therm_eff;
03268 break;
03269 case STOGNIENKO95:
03270 case BRUGGEMAN35:
03271
03272 gd->work = MIN2(gd->work,gdArr[i].work);
03273 gd->bandgap = MIN2(gd->bandgap,gdArr[i].bandgap);
03274 gd->therm_eff += frac2[i]*gdArr[i].therm_eff;
03275 break;
03276 default:
03277 fprintf( ioQQQ, " Insanity in mie_read_mix\n" );
03278 ShowMe();
03279 puts( "[Stop in mie_read_mix]" );
03280 cdEXIT(EXIT_FAILURE);
03281 }
03282 gd->matType = gdArr[i].matType;
03283 gd->subl_temp = MIN2(gd->subl_temp,gdArr[i].subl_temp);
03284 }
03285 }
03286
03287 if( gd->rho <= 0. )
03288 {
03289 fprintf( ioQQQ, " Illegal value for the density: %.3e\n", gd->rho );
03290 puts( "[Stop in mie_read_mix]" );
03291 cdEXIT(EXIT_FAILURE);
03292 }
03293 if( gd->mol_weight <= 0. )
03294 {
03295 fprintf( ioQQQ, " Illegal value for the molecular weight: %.3e\n", gd->mol_weight );
03296 puts( "[Stop in mie_read_mix]" );
03297 cdEXIT(EXIT_FAILURE);
03298 }
03299 if( maxIndex <= 0. )
03300 {
03301 fprintf( ioQQQ, " No atoms were found in the grain molecule\n" );
03302 puts( "[Stop in mie_read_mix]" );
03303 cdEXIT(EXIT_FAILURE);
03304 }
03305
03306
03307 ASSERT( wavLo > 0. && wavHi < DBL_MAX && wavLo < wavHi );
03308 ASSERT( gd->abun > 0. && gd->abun < DBL_MAX );
03309 ASSERT( gd->work > 0. && gd->work < DBL_MAX );
03310 ASSERT( gd->bandgap >= 0. && gd->bandgap < gd->work );
03311 ASSERT( gd->therm_eff > 0. && gd->therm_eff <= 1. );
03312 ASSERT( gd->subl_temp > 0. && gd->subl_temp < DBL_MAX );
03313 ASSERT( minIndex > 0. && minIndex < DBL_MAX );
03314
03315
03316 wavLo *= 1. + 10.*DBL_EPSILON;
03317 wavHi *= 1. - 10.*DBL_EPSILON;
03318
03319
03320 nAtoms = 0.;
03321 for( nelem=0; nelem < LIMELM; nelem++ )
03322 {
03323 gd->elmAbun[nelem] /= minIndex;
03324 nAtoms += gd->elmAbun[nelem];
03325 }
03326 ASSERT( nAtoms > 0. );
03327 gd->abun *= minIndex;
03328 gd->mol_weight /= minIndex;
03329
03330 gd->atom_weight = gd->mol_weight/nAtoms;
03331
03332 mie_write_form(gd->elmAbun,chWord);
03333 fprintf( ioQQQ, "\n The chemical formula of the new grain molecule is: %s\n", chWord );
03334 fprintf( ioQQQ, " The abundance wrt H at maximum depletion of this molecule is: %.3e\n",
03335 gd->abun );
03336 fprintf( ioQQQ, " The abundance wrt H at standard depletion of this molecule is: %.3e\n",
03337 gd->abun*gd->depl );
03338
03339
03340 for( nelem=0; nelem < LIMELM; nelem++ )
03341 {
03342 gd->elmAbun[nelem] *= gd->abun*gd->depl;
03343 }
03344
03345 delta = (double*)MALLOC(sizeof(double)*(unsigned)sumAxes);
03346 frdelta = (double*)MALLOC(sizeof(double)*(unsigned)sumAxes);
03347 eps = (complex<double>*)MALLOC(sizeof(complex<double>)*(unsigned)sumAxes);
03348
03349 l = 0;
03350 for( i=0; i < nMaterial; i++ )
03351 {
03352 for( k=0; k < gdArr[i].nAxes; k++ )
03353 {
03354 frdelta[l] = gdArr[i].wt[k]*frac[i];
03355 delta[l] = ( l == 0 ) ? frdelta[l] : delta[l-1] + frdelta[l];
03356 ++l;
03357 }
03358 }
03359 ASSERT( l == sumAxes && fabs(delta[l-1]-1.) < 10.*DBL_EPSILON );
03360
03361
03362
03363 gd->wavlen[0] = (double*)MALLOC(sizeof(double)*(unsigned)gd->ndata[0]);
03364 gd->n[0] = (complex<double>*)MALLOC(sizeof(complex<double>)*(unsigned)gd->ndata[0]);
03365 gd->nr1[0] = (double*)MALLOC(sizeof(double)*(unsigned)gd->ndata[0]);
03366
03367 wavStep = log(wavHi/wavLo)/(double)(gd->ndata[0]-1);
03368
03369 switch( EMTtype )
03370 {
03371 case FARAFONOV00:
03372
03373
03374
03375
03376
03377
03378 for( j=0; j < gd->ndata[0]; j++ )
03379 {
03380 double nre,nim;
03381 complex<double> a1,a2,a1c,a2c,a11,a12,a21,a22,ratio;
03382
03383 gd->wavlen[0][j] = wavLo*exp((double)j*wavStep);
03384
03385 init_eps(gd->wavlen[0][j],nMaterial,gdArr,eps);
03386
03387 ratio = eps[0]/eps[1];
03388
03389 a1 = (ratio+2.)/3.;
03390 a2 = (1.-ratio)*delta[0];
03391
03392 for( l=1; l < sumAxes-1; l++ )
03393 {
03394 ratio = eps[l]/eps[l+1];
03395
03396 a1c = a1;
03397 a2c = a2;
03398 a11 = (ratio+2.)/3.;
03399 a12 = (2.-2.*ratio)/(9.*delta[l]);
03400 a21 = (1.-ratio)*delta[l];
03401 a22 = (2.*ratio+1.)/3.;
03402
03403 a1 = a11*a1c + a12*a2c;
03404 a2 = a21*a1c + a22*a2c;
03405 }
03406
03407 a1c = a1;
03408 a2c = a2;
03409 a11 = 1.;
03410 a12 = 1./3.;
03411 a21 = eps[sumAxes-1];
03412 a22 = -2./3.*eps[sumAxes-1];
03413
03414 a1 = a11*a1c + a12*a2c;
03415 a2 = a21*a1c + a22*a2c;
03416
03417 ratio = a2/a1;
03418 dftori(&nre,&nim,ratio.real(),ratio.imag());
03419
03420 gd->n[0][j] = complex<double>(nre,nim);
03421 gd->nr1[0][j] = nre-1.;
03422 }
03423 break;
03424 case STOGNIENKO95:
03425 case BRUGGEMAN35:
03426 for( j=0; j < gd->ndata[0]; j++ )
03427 {
03428 const double EPS_TOLER = 10.*DBL_EPSILON;
03429 double nre,nim;
03430 complex<double> eps0;
03431
03432 gd->wavlen[0][j] = wavLo*exp((double)j*wavStep);
03433
03434 init_eps(gd->wavlen[0][j],nMaterial,gdArr,eps);
03435
03436
03437 if( j == 0 )
03438 {
03439
03440 eps0 = 0.;
03441 for( l=0; l < sumAxes; l++ )
03442 eps0 += frdelta[l]*eps[l];
03443 }
03444 else
03445 {
03446
03447 eps0 = eps_eff;
03448 }
03449
03450 if( EMTtype == STOGNIENKO95 )
03451
03452
03453 eps_eff = cnewton( Stognienko, frdelta, eps, sumAxes, eps0, EPS_TOLER );
03454 else if( EMTtype == BRUGGEMAN35 )
03455
03456
03457 eps_eff = cnewton( Bruggeman, frdelta, eps, sumAxes, eps0, EPS_TOLER );
03458 else
03459 {
03460 fprintf( ioQQQ, " Insanity in mie_read_mix\n" );
03461 ShowMe();
03462 puts( "[Stop in mie_read_mix]" );
03463 cdEXIT(EXIT_FAILURE);
03464 }
03465
03466 dftori(&nre,&nim,eps_eff.real(),eps_eff.imag());
03467
03468 gd->n[0][j] = complex<double>(nre,nim);
03469 gd->nr1[0][j] = nre-1.;
03470 }
03471 break;
03472 default:
03473 fprintf( ioQQQ, " Insanity in mie_read_mix\n" );
03474 ShowMe();
03475 puts( "[Stop in mie_read_mix]" );
03476 cdEXIT(EXIT_FAILURE);
03477 }
03478
03479
03480 for( i=0; i < nMaterial; i++ )
03481 {
03482 for( j=0; j < gdArr[i].nOpcCols; j++ )
03483 {
03484 if( gdArr[i].opcData[j] != NULL )
03485 free(gdArr[i].opcData[j]);
03486 }
03487 if( gdArr[i].opcAnu != NULL )
03488 free(gdArr[i].opcAnu);
03489 for( j=0; j < gdArr[i].nAxes; j++ )
03490 {
03491 if( gdArr[i].nr1[j] != NULL )
03492 free(gdArr[i].nr1[j]);
03493 if( gdArr[i].n[j] != NULL )
03494 free(gdArr[i].n[j]);
03495 if( gdArr[i].wavlen[j] != NULL )
03496 free(gdArr[i].wavlen[j]);
03497 }
03498 }
03499
03500 free(eps);
03501 free(frdelta);
03502 free(delta);
03503 free(gdArr);
03504 free(frac2);
03505 free(frac);
03506
03507 DEBUG_EXIT( "mie_read_mix()" );
03508 return;
03509 }
03510
03511
03512
03513 static void init_eps(double wavlen,
03514 long nMaterial,
03515 grain_data gdArr[],
03516 complex<double> eps[])
03517 {
03518 long i,
03519 k;
03520
03521 long l = 0;
03522
03523 DEBUG_ENTRY( "init_eps()" );
03524 for( i=0; i < nMaterial; i++ )
03525 {
03526 for( k=0; k < gdArr[i].nAxes; k++ )
03527 {
03528 bool lgErr;
03529 long ind;
03530 double eps1,eps2,frc,nim,nre;
03531
03532 find_arr(wavlen,gdArr[i].wavlen[k],gdArr[i].ndata[k],&ind,&lgErr);
03533 ASSERT( !lgErr );
03534 frc = (wavlen-gdArr[i].wavlen[k][ind])/(gdArr[i].wavlen[k][ind+1]-gdArr[i].wavlen[k][ind]);
03535 ASSERT( frc > 0.-10.*DBL_EPSILON && frc < 1.+10.*DBL_EPSILON );
03536 nre = (1.-frc)*gdArr[i].n[k][ind].real() + frc*gdArr[i].n[k][ind+1].real();
03537 ASSERT( nre > 0. );
03538 nim = (1.-frc)*gdArr[i].n[k][ind].imag() + frc*gdArr[i].n[k][ind+1].imag();
03539 ASSERT( nim >= 0. );
03540 ritodf(nre,nim,&eps1,&eps2);
03541 eps[l++] = complex<double>(eps1,eps2);
03542 }
03543 }
03544
03545 DEBUG_EXIT( "init_eps()" );
03546 return;
03547 }
03548
03549
03550
03551
03552
03553
03554
03555
03556
03557
03558
03559 static complex<double> cnewton(
03560 void(*fun)(complex<double>,double[],complex<double>[],long,complex<double>*,double*,double*),
03561 double frdelta[],
03562 complex<double> eps[],
03563 long sumAxes,
03564 complex<double> x0,
03565 double tol)
03566 {
03567 int i;
03568
03569 const int LOOP_MAX = 100;
03570 const double TINY = 1.e-12;
03571
03572 DEBUG_ENTRY( "cnewton()" );
03573 for( i=0; i < LOOP_MAX; i++ )
03574 {
03575 complex<double> x1,y;
03576 double dudx,dudy,norm2;
03577
03578 (*fun)(x0,frdelta,eps,sumAxes,&y,&dudx,&dudy);
03579
03580 norm2 = POW2(dudx) + POW2(dudy);
03581
03582 if( norm2 < TINY*norm(y) )
03583 {
03584 fprintf( ioQQQ, " cnewton - zero divide error\n" );
03585 ShowMe();
03586 puts( "[Stop in cnewton]" );
03587 cdEXIT(EXIT_FAILURE);
03588 }
03589 x1 = x0 - complex<double>( y.real()*dudx-y.imag()*dudy, y.imag()*dudx+y.real()*dudy )/norm2;
03590
03591
03592 if( fabs(x0.real()/x1.real()-1.) + fabs(x0.imag()/x1.imag()-1.) < tol )
03593 {
03594 DEBUG_EXIT( "cnewton()" );
03595 return x1;
03596 }
03597
03598 x0 = x1;
03599 }
03600
03601 fprintf( ioQQQ, " cnewton did not converge\n" );
03602 ShowMe();
03603 puts( "[Stop in cnewton]" );
03604 cdEXIT(EXIT_FAILURE);
03605 }
03606
03607
03608
03609
03610
03611 static void Stognienko(complex<double> x,
03612 double frdelta[],
03613 complex<double> eps[],
03614 long sumAxes,
03615 complex<double> *f,
03616 double *dudx,
03617 double *dudy)
03618 {
03619 long i,
03620 l;
03621
03622 static const double L[4] = { 0., 1./2., 1., 1./3. };
03623 static const double fl[4] = { 5./9., 2./9., 2./9., 1. };
03624
03625 DEBUG_ENTRY( "Stognienko()" );
03626 *f = complex<double>(0.,0.);
03627 *dudx = 0.;
03628 *dudy = 0.;
03629 for( l=0; l < sumAxes; l++ )
03630 {
03631 complex<double> hlp = eps[l] - x;
03632 double h1 = eps[l].imag()*x.real() - eps[l].real()*x.imag();
03633
03634 for( i=0; i < 4; i++ )
03635 {
03636 double f1 = fl[i]*frdelta[l];
03637 double xx = ( i < 3 ) ? sin(PI*frdelta[l]) : cos(PI*frdelta[l]);
03638 complex<double> f2 = f1*xx*xx;
03639 complex<double> one = x + hlp*L[i];
03640 complex<double> two = f2*hlp/one;
03641 double h2 = norm(one);
03642 *f += two;
03643 *dudx -= f2.real()*(eps[l].real()*h2 + h1*2.*one.imag()*(1.-L[i]))/POW2(h2);
03644 *dudy -= f2.real()*(eps[l].imag()*h2 - h1*2.*one.real()*(1.-L[i]))/POW2(h2);
03645 }
03646 }
03647
03648 DEBUG_EXIT( "Stognienko()" );
03649 return;
03650 }
03651
03652
03653
03654
03655 static void Bruggeman(complex<double> x,
03656 double frdelta[],
03657 complex<double> eps[],
03658 long sumAxes,
03659 complex<double> *f,
03660 double *dudx,
03661 double *dudy)
03662 {
03663 long l;
03664
03665 static const double L = 1./3.;
03666
03667 DEBUG_ENTRY( "Bruggeman()" );
03668 *f = complex<double>(0.,0.);
03669 *dudx = 0.;
03670 *dudy = 0.;
03671 for( l=0; l < sumAxes; l++ )
03672 {
03673 complex<double> hlp = eps[l] - x;
03674 double h1 = eps[l].imag()*x.real() - eps[l].real()*x.imag();
03675 complex<double> f2 = frdelta[l];
03676 complex<double> one = x + hlp*L;
03677 complex<double> two = f2*hlp/one;
03678 double h2 = norm(one);
03679 *f += two;
03680 *dudx -= f2.real()*(eps[l].real()*h2 + h1*2.*one.imag()*(1.-L))/POW2(h2);
03681 *dudy -= f2.real()*(eps[l].imag()*h2 - h1*2.*one.real()*(1.-L))/POW2(h2);
03682 }
03683
03684 DEBUG_EXIT( "Bruggeman()" );
03685 return;
03686 }
03687
03688
03689
03690 static void mie_read_szd( const char *chFile,
03691 sd_data *sd)
03692 {
03693 bool lgTryOverride = false;
03694 long int dl,
03695 i;
03696 double mul = 0.,
03697 ref_neg = DBL_MAX,
03698 ref_pos = DBL_MAX,
03699 step_neg = DBL_MAX,
03700 step_pos = DBL_MAX;
03701 char chLine[FILENAME_PATH_LENGTH_2],
03702 chWord[WORDLEN];
03703 FILE *io2;
03704
03705
03706
03707
03708
03709
03710 const double FRAC_CUTOFF = 1.e-4;
03711 const double MUL_CO1 = -log(FRAC_CUTOFF);
03712 const double MUL_CO2 = sqrt(MUL_CO1);
03713 const double MUL_CO3 = pow(MUL_CO1,1./3.);
03714 const double MUL_LND = sqrt(-2.*log(FRAC_CUTOFF));
03715 const double MUL_NRM = MUL_LND;
03716
03717 DEBUG_ENTRY( "mie_read_szd()" );
03718 if( (io2 = fopen(chFile,"r")) == NULL )
03719 {
03720 fprintf( ioQQQ, " Could not open %s for reading\n",chFile);
03721 puts( "[Stop in mie_read_szd]" );
03722 cdEXIT(EXIT_FAILURE);
03723 }
03724
03725 dl = 0;
03726
03727
03728 mie_next_data(chFile,io2,chLine,&dl);
03729 mie_read_long(chFile,chLine,&sd->magic,true,dl);
03730 if( sd->magic != MAGIC_SZD )
03731 {
03732 fprintf( ioQQQ, " Size distribution file %s has obsolete magic number\n",chFile );
03733 fprintf( ioQQQ, " I found %ld, but expected %ld on line #%ld\n",sd->magic,MAGIC_SZD,dl );
03734 fprintf( ioQQQ, " Please replace this file with an up to date version\n" );
03735 puts( "[Stop in mie_read_szd]" );
03736 cdEXIT(EXIT_FAILURE);
03737 }
03738
03739
03740 mie_next_data(chFile,io2,chLine,&dl);
03741 mie_read_word(chLine,chWord,WORDLEN,true);
03742
03743 if( nMatch( "SSIZ", chWord ) )
03744 {
03745 sd->sdCase = SD_SINGLE_SIZE;
03746 }
03747 else if( nMatch( "POWE", chWord ) )
03748 {
03749 sd->sdCase = SD_POWERLAW;
03750 }
03751 else if( nMatch( "EXP1", chWord ) )
03752 {
03753 sd->sdCase = SD_EXP_CUTOFF1;
03754 sd->a[ipAlpha] = 1.;
03755 mul = MUL_CO1;
03756 }
03757 else if( nMatch( "EXP2", chWord ) )
03758 {
03759 sd->sdCase = SD_EXP_CUTOFF2;
03760 sd->a[ipAlpha] = 2.;
03761 mul = MUL_CO2;
03762 }
03763 else if( nMatch( "EXP3", chWord ) )
03764 {
03765 sd->sdCase = SD_EXP_CUTOFF3;
03766 sd->a[ipAlpha] = 3.;
03767 mul = MUL_CO3;
03768 }
03769 else if( nMatch( "LOGN", chWord ) )
03770 {
03771 sd->sdCase = SD_LOG_NORMAL;
03772 mul = MUL_LND;
03773 }
03774
03775 else if( nMatch( "NORM", chWord ) )
03776 {
03777 sd->sdCase = SD_LIN_NORMAL;
03778 mul = MUL_NRM;
03779 }
03780 else if( nMatch( "TABL", chWord ) )
03781 {
03782 sd->sdCase = SD_TABLE;
03783 }
03784 else
03785 {
03786 sd->sdCase = SD_ILLEGAL;
03787 }
03788
03789 switch( sd->sdCase )
03790 {
03791 case SD_SINGLE_SIZE:
03792
03793 mie_next_data(chFile,io2,chLine,&dl);
03794 mie_read_double(chFile,chLine,&sd->a[ipSize],true,dl);
03795 if( sd->a[ipSize] < SMALLEST_GRAIN || sd->a[ipSize] > LARGEST_GRAIN )
03796 {
03797 fprintf( ioQQQ, " Illegal value for grain size\n" );
03798 fprintf( ioQQQ, " Grain sizes should be between %.5f and %.0f micron\n",
03799 SMALLEST_GRAIN, LARGEST_GRAIN );
03800 fprintf( ioQQQ, " Line #%ld: %s\n",dl,chLine);
03801 puts( "[Stop in mie_read_szd]" );
03802 cdEXIT(EXIT_FAILURE);
03803 }
03804 break;
03805 case SD_POWERLAW:
03806
03807 mie_next_data(chFile,io2,chLine,&dl);
03808 mie_read_double(chFile,chLine,&sd->a[ipBLo],true,dl);
03809 if( sd->a[ipBLo] < SMALLEST_GRAIN || sd->a[ipBLo] > LARGEST_GRAIN )
03810 {
03811 fprintf( ioQQQ, " Illegal value for grain size (lower limit)\n" );
03812 fprintf( ioQQQ, " Grain sizes should be between %.5f and %.0f micron\n",
03813 SMALLEST_GRAIN, LARGEST_GRAIN );
03814 fprintf( ioQQQ, " Line #%ld: %s\n",dl,chLine);
03815 puts( "[Stop in mie_read_szd]" );
03816 cdEXIT(EXIT_FAILURE);
03817 }
03818
03819
03820 mie_next_data(chFile,io2,chLine,&dl);
03821 mie_read_double(chFile,chLine,&sd->a[ipBHi],true,dl);
03822 if( sd->a[ipBHi] < SMALLEST_GRAIN || sd->a[ipBHi] > LARGEST_GRAIN ||
03823 sd->a[ipBHi] <= sd->a[ipBLo] )
03824 {
03825 fprintf( ioQQQ, " Illegal value for grain size (upper limit)\n" );
03826 fprintf( ioQQQ, " Grain sizes should be between %.5f and %.0f micron\n",
03827 SMALLEST_GRAIN, LARGEST_GRAIN );
03828 fprintf( ioQQQ, " and upper limit should be greater than lower limit\n" );
03829 fprintf( ioQQQ, " Line #%ld: %s\n",dl,chLine);
03830 puts( "[Stop in mie_read_szd]" );
03831 cdEXIT(EXIT_FAILURE);
03832 }
03833
03834
03835 mie_next_data(chFile,io2,chLine,&dl);
03836 if( sscanf( chLine, "%lf", &sd->a[ipExp] ) != 1 )
03837 {
03838 fprintf( ioQQQ, " Syntax error in %s\n",chFile);
03839 fprintf( ioQQQ, " Line #%ld: %s\n",dl,chLine);
03840 puts( "[Stop in mie_read_szd]" );
03841 cdEXIT(EXIT_FAILURE);
03842 }
03843
03844 sd->a[ipBeta] = 0.;
03845 sd->a[ipSLo] = 0.;
03846 sd->a[ipSHi] = 0.;
03847
03848 sd->lim[ipBLo] = sd->a[ipBLo];
03849 sd->lim[ipBHi] = sd->a[ipBHi];
03850 break;
03851 case SD_EXP_CUTOFF1:
03852 case SD_EXP_CUTOFF2:
03853 case SD_EXP_CUTOFF3:
03854
03855
03856
03857 mie_next_data(chFile,io2,chLine,&dl);
03858 mie_read_double(chFile,chLine,&sd->a[ipBLo],true,dl);
03859
03860
03861 mie_next_data(chFile,io2,chLine,&dl);
03862 mie_read_double(chFile,chLine,&sd->a[ipBHi],true,dl);
03863
03864
03865 mie_next_data(chFile,io2,chLine,&dl);
03866 if( sscanf( chLine, "%lf", &sd->a[ipExp] ) != 1 )
03867 {
03868 fprintf( ioQQQ, " Syntax error in %s\n",chFile);
03869 fprintf( ioQQQ, " Line #%ld: %s\n",dl,chLine);
03870 puts( "[Stop in mie_read_szd]" );
03871 cdEXIT(EXIT_FAILURE);
03872 }
03873
03874
03875 mie_next_data(chFile,io2,chLine,&dl);
03876 if( sscanf( chLine, "%lf", &sd->a[ipBeta] ) != 1 )
03877 {
03878 fprintf( ioQQQ, " Syntax error in %s\n",chFile);
03879 fprintf( ioQQQ, " Line #%ld: %s\n",dl,chLine);
03880 puts( "[Stop in mie_read_szd]" );
03881 cdEXIT(EXIT_FAILURE);
03882 }
03883
03884
03885 mie_next_data(chFile,io2,chLine,&dl);
03886 mie_read_double(chFile,chLine,&sd->a[ipSLo],false,dl);
03887
03888
03889 mie_next_data(chFile,io2,chLine,&dl);
03890 mie_read_double(chFile,chLine,&sd->a[ipSHi],false,dl);
03891
03892 ref_neg = sd->a[ipBLo];
03893 step_neg = -mul*sd->a[ipSLo];
03894 ref_pos = sd->a[ipBHi];
03895 step_pos = mul*sd->a[ipSHi];
03896 lgTryOverride = true;
03897 break;
03898 case SD_LOG_NORMAL:
03899
03900 mie_next_data(chFile,io2,chLine,&dl);
03901 mie_read_double(chFile,chLine,&sd->a[ipGCen],true,dl);
03902
03903
03904 mie_next_data(chFile,io2,chLine,&dl);
03905 mie_read_double(chFile,chLine,&sd->a[ipGSig],true,dl);
03906
03907
03908 ref_neg = ref_pos = sd->a[ipGCen]*exp(3.*POW2(sd->a[ipGSig]));
03909 step_neg = sd->a[ipGCen]*(exp(-mul*sd->a[ipGSig]) - 1.);
03910 step_pos = sd->a[ipGCen]*(exp(mul*sd->a[ipGSig]) - 1.);
03911 lgTryOverride = true;
03912 break;
03913 case SD_LIN_NORMAL:
03914
03915 mie_next_data(chFile,io2,chLine,&dl);
03916 mie_read_double(chFile,chLine,&sd->a[ipGCen],true,dl);
03917
03918
03919 mie_next_data(chFile,io2,chLine,&dl);
03920 mie_read_double(chFile,chLine,&sd->a[ipGSig],true,dl);
03921
03922
03923 ref_neg = ref_pos = (sd->a[ipGCen] + sqrt(POW2(sd->a[ipGCen]) + 12.*POW2(sd->a[ipGSig])))/2.;
03924 step_neg = -mul*sd->a[ipGSig];
03925 step_pos = mul*sd->a[ipGSig];
03926 lgTryOverride = true;
03927 break;
03928 case SD_TABLE:
03929
03930 mie_next_data(chFile,io2,chLine,&dl);
03931 mie_read_double(chFile,chLine,&sd->a[ipBLo],true,dl);
03932 if( sd->a[ipBLo] < SMALLEST_GRAIN || sd->a[ipBLo] > LARGEST_GRAIN )
03933 {
03934 fprintf( ioQQQ, " Illegal value for grain size (lower limit)\n" );
03935 fprintf( ioQQQ, " Grain sizes should be between %.5f and %.0f micron\n",
03936 SMALLEST_GRAIN, LARGEST_GRAIN );
03937 fprintf( ioQQQ, " Line #%ld: %s\n",dl,chLine);
03938 puts( "[Stop in mie_read_szd]" );
03939 cdEXIT(EXIT_FAILURE);
03940 }
03941
03942
03943 mie_next_data(chFile,io2,chLine,&dl);
03944 mie_read_double(chFile,chLine,&sd->a[ipBHi],true,dl);
03945 if( sd->a[ipBHi] < SMALLEST_GRAIN || sd->a[ipBHi] > LARGEST_GRAIN ||
03946 sd->a[ipBHi] <= sd->a[ipBLo] )
03947 {
03948 fprintf( ioQQQ, " Illegal value for grain size (upper limit)\n" );
03949 fprintf( ioQQQ, " Grain sizes should be between %.5f and %.0f micron\n",
03950 SMALLEST_GRAIN, LARGEST_GRAIN );
03951 fprintf( ioQQQ, " and upper limit should be greater than lower limit\n" );
03952 fprintf( ioQQQ, " Line #%ld: %s\n",dl,chLine);
03953 puts( "[Stop in mie_read_szd]" );
03954 cdEXIT(EXIT_FAILURE);
03955 }
03956
03957
03958 mie_next_data(chFile,io2,chLine,&dl);
03959 mie_read_long(chFile,chLine,&sd->npts,true,dl);
03960 if( sd->npts < 2 )
03961 {
03962 fprintf( ioQQQ, " Illegal value for no. of points in table\n" );
03963 fprintf( ioQQQ, " Line #%ld: %s\n",dl,chLine);
03964 puts( "[Stop in mie_read_szd]" );
03965 cdEXIT(EXIT_FAILURE);
03966 }
03967
03968
03969 sd->ln_a = (double *)MALLOC(sizeof(double)*(unsigned)sd->npts);
03970 sd->ln_a4dNda = (double *)MALLOC(sizeof(double)*(unsigned)sd->npts);
03971 if( sd->ln_a == NULL || sd->ln_a4dNda == NULL )
03972 {
03973 fprintf( ioQQQ, " Could not MALLOC arrays\n" );
03974 puts( "[Stop in mie_read_szd]" );
03975 cdEXIT(EXIT_FAILURE);
03976 }
03977
03978
03979 for( i=0; i < sd->npts; ++i )
03980 {
03981 double help1, help2;
03982
03983 mie_next_data(chFile,io2,chLine,&dl);
03984
03985 if( sscanf( chLine, "%le %le", &help1, &help2 ) != 2 )
03986 {
03987 fprintf( ioQQQ, " Syntax error in %s\n",chFile);
03988 fprintf( ioQQQ, " Line #%ld: %s\n",dl,chLine);
03989 puts( "[Stop in mie_read_szd]" );
03990 cdEXIT(EXIT_FAILURE);
03991 }
03992
03993 if( help1 <= 0. || help2 <= 0. )
03994 {
03995 fprintf( ioQQQ, " Reading table failed on line #%ld of %s\n",dl,chFile );
03996 fprintf( ioQQQ, " Illegal data value %.6e or %.6e\n", help1, help2 );
03997 puts( "[Stop in mie_read_szd]" );
03998 cdEXIT(EXIT_FAILURE);
03999 }
04000
04001 sd->ln_a[i] = log(help1);
04002 sd->ln_a4dNda[i] = log(help2);
04003
04004 if( i > 0 && sd->ln_a[i] <= sd->ln_a[i-1] )
04005 {
04006 fprintf( ioQQQ, " Reading table failed on line #%ld of %s\n",dl,chFile );
04007 fprintf( ioQQQ, " Grain radii should be monotonically increasing\n" );
04008 puts( "[Stop in mie_read_szd]" );
04009 cdEXIT(EXIT_FAILURE);
04010 }
04011 }
04012
04013 sd->lim[ipBLo] = sd->a[ipBLo];
04014 sd->lim[ipBHi] = sd->a[ipBHi];
04015 break;
04016 case SD_ILLEGAL:
04017 default:
04018 fprintf( ioQQQ, " unimplemented case for grain size distribution in file %s\n" , chFile );
04019 fprintf( ioQQQ, " Line #%ld: value %s\n",dl,chWord);
04020 puts( "[Stop in mie_read_szd]" );
04021 cdEXIT(EXIT_FAILURE);
04022 }
04023
04024
04025
04026
04027
04028
04029
04030 if( lgTryOverride )
04031 {
04032 double help;
04033
04034 mie_next_data(chFile,io2,chLine,&dl);
04035 mie_read_double(chFile,chLine,&help,false,dl);
04036 sd->lim[ipBLo] = ( help <= 0. ) ? search_limit(ref_neg,step_neg,FRAC_CUTOFF,*sd) : help;
04037
04038 mie_next_data(chFile,io2,chLine,&dl);
04039 mie_read_double(chFile,chLine,&help,false,dl);
04040 sd->lim[ipBHi] = ( help <= 0. ) ? search_limit(ref_pos,step_pos,FRAC_CUTOFF,*sd) : help;
04041
04042 if( sd->lim[ipBLo] < SMALLEST_GRAIN || sd->lim[ipBHi] > LARGEST_GRAIN ||
04043 sd->lim[ipBHi] <= sd->lim[ipBLo] )
04044 {
04045 fprintf( ioQQQ, " Illegal size limits: lower %.5f and/or upper %.5f\n",
04046 sd->lim[ipBLo], sd->lim[ipBHi] );
04047 fprintf( ioQQQ, " Grain sizes should be between %.5f and %.0f micron\n",
04048 SMALLEST_GRAIN, LARGEST_GRAIN );
04049 fprintf( ioQQQ, " and upper limit should be greater than lower limit\n" );
04050 fprintf( ioQQQ, " Please alter the size distribution file\n" );
04051 puts( "[Stop in mie_read_szd]" );
04052 cdEXIT(EXIT_FAILURE);
04053 }
04054 }
04055
04056 fclose(io2);
04057
04058 DEBUG_EXIT( "mie_read_szd()" );
04059 return;
04060 }
04061
04062
04063 static void mie_read_long( const char *chFile,
04064 const char chLine[],
04065 long int *data,
04066 bool lgZeroIllegal,
04067 long int dl)
04068 {
04069 DEBUG_ENTRY( "mie_read_long()" );
04070 if( sscanf( chLine, "%ld", data ) != 1 )
04071 {
04072 fprintf( ioQQQ, " Syntax error in %s\n",chFile);
04073 fprintf( ioQQQ, " Line #%ld: %s\n",dl,chLine);
04074 puts( "[Stop in mie_read_long]" );
04075 cdEXIT(EXIT_FAILURE);
04076 }
04077 if( *data < 0 || (*data == 0 && lgZeroIllegal) )
04078 {
04079 fprintf( ioQQQ, " Illegal data value in %s\n",chFile);
04080 fprintf( ioQQQ, " Line #%ld: %ld\n",dl,*data);
04081 puts( "[Stop in mie_read_long]" );
04082 cdEXIT(EXIT_FAILURE);
04083 }
04084 DEBUG_EXIT( "mie_read_long()" );
04085 return;
04086 }
04087
04088
04089 static void mie_read_float( const char *chFile,
04090 const char chLine[],
04091 float *data,
04092 bool lgZeroIllegal,
04093 long int dl)
04094 {
04095 DEBUG_ENTRY( "mie_read_float()" );
04096 if( sscanf( chLine, "%f", data ) != 1 )
04097 {
04098 fprintf( ioQQQ, " Syntax error in %s\n",chFile);
04099 fprintf( ioQQQ, " Line #%ld: %s\n",dl,chLine);
04100 puts( "[Stop in mie_read_float]" );
04101 cdEXIT(EXIT_FAILURE);
04102 }
04103 if( *data < 0. || (*data == 0. && lgZeroIllegal) )
04104 {
04105 fprintf( ioQQQ, " Illegal data value in %s\n",chFile);
04106 fprintf( ioQQQ, " Line #%ld: %14.6e\n",dl,*data);
04107 puts( "[Stop in mie_read_float]" );
04108 cdEXIT(EXIT_FAILURE);
04109 }
04110 DEBUG_EXIT( "mie_read_float()" );
04111 return;
04112 }
04113
04114
04115 static void mie_read_double( const char *chFile,
04116 const char chLine[],
04117 double *data,
04118 bool lgZeroIllegal,
04119 long int dl)
04120 {
04121 DEBUG_ENTRY( "mie_read_double()" );
04122 if( sscanf( chLine, "%lf", data ) != 1 )
04123 {
04124 fprintf( ioQQQ, " Syntax error in %s\n",chFile);
04125 fprintf( ioQQQ, " Line #%ld: %s\n",dl,chLine);
04126 puts( "[Stop in mie_read_double]" );
04127 cdEXIT(EXIT_FAILURE);
04128 }
04129 if( *data < 0. || (*data == 0. && lgZeroIllegal) )
04130 {
04131 fprintf( ioQQQ, " Illegal data value in %s\n",chFile);
04132 fprintf( ioQQQ, " Line #%ld: %14.6e\n",dl,*data);
04133 puts( "[Stop in mie_read_double]" );
04134 cdEXIT(EXIT_FAILURE);
04135 }
04136 DEBUG_EXIT( "mie_read_double()" );
04137 return;
04138 }
04139
04140
04141 static void mie_read_form( const char chWord[],
04142 double elmAbun[],
04143 double *no_atoms,
04144 double *mol_weight)
04145 {
04146 long int nelem,
04147 len;
04148 char chElmName[3];
04149 double frac;
04150
04151 DEBUG_ENTRY( "mie_read_form()" );
04152
04153 *no_atoms = 0.;
04154 *mol_weight = 0.;
04155 string strWord( chWord );
04156 for( nelem=0; nelem < LIMELM; nelem++ )
04157 {
04158 frac = 0.;
04159 strcpy(chElmName,elementnames.chElementSym[nelem]);
04160 if( chElmName[1] == ' ' )
04161 chElmName[1] = '\0';
04162 string::size_type ptr = strWord.find( chElmName );
04163 if( ptr != string::npos )
04164 {
04165 len = (long)strlen(chElmName);
04166
04167 if( !islower((unsigned char)chWord[ptr+len]) )
04168 {
04169 if( isdigit((unsigned char)chWord[ptr+len]) )
04170 {
04171 sscanf(chWord+ptr+len,"%lf",&frac);
04172 }
04173 else
04174 {
04175 frac = 1.;
04176 }
04177 }
04178 }
04179 elmAbun[nelem] = frac;
04180
04181 if( nelem != ipHYDROGEN )
04182 *no_atoms += frac;
04183 *mol_weight += frac*dense.AtomicWeight[nelem];
04184 }
04185
04186 if( *no_atoms == 0. )
04187 *no_atoms = 1.;
04188
04189 DEBUG_EXIT( "mie_read_form()" );
04190 return;
04191 }
04192
04193
04194 static void mie_write_form( const double elmAbun[],
04195 char chWord[])
04196 {
04197 long int len,
04198 nelem;
04199
04200 DEBUG_ENTRY( "mie_write_form()" );
04201
04202 len=0;
04203 for( nelem=0; nelem < LIMELM; nelem++ )
04204 {
04205 if( elmAbun[nelem] > 0. )
04206 {
04207 char chElmName[3];
04208 long index100 = nint(100.*elmAbun[nelem]);
04209
04210 strcpy(chElmName,elementnames.chElementSym[nelem]);
04211 if( chElmName[1] == ' ' )
04212 chElmName[1] = '\0';
04213
04214 if( index100 == 100 )
04215 sprintf( &chWord[len], "%s", chElmName );
04216 else if( index100%100 == 0 )
04217 sprintf( &chWord[len], "%s%ld", chElmName, index100/100 );
04218 else
04219 {
04220 double xIndex = (double)index100/100.;
04221 sprintf( &chWord[len], "%s%.2f", chElmName, xIndex );
04222 }
04223 len = (long)strlen( chWord );
04224 ASSERT( len < FILENAME_PATH_LENGTH_2-25 );
04225 }
04226 }
04227
04228 DEBUG_EXIT( "mie_write_form()" );
04229 return;
04230 }
04231
04232
04233 static void mie_read_word( const char chLine[],
04234 char chWord[],
04235 long n,
04236 int toUpper)
04237 {
04238 long ip = 0,
04239 op = 0;
04240
04241 DEBUG_ENTRY( "mie_read_word()" );
04242
04243
04244 while( chLine[ip] == ' ' || chLine[ip] == '"' )
04245 ip++;
04246
04247 while( op < n-1 && chLine[ip] != ' ' && chLine[ip] != '"' )
04248 if( toUpper )
04249 chWord[op++] = (char)toupper(chLine[ip++]);
04250 else
04251 chWord[op++] = chLine[ip++];
04252
04253 chWord[op] = '\0';
04254
04255 DEBUG_EXIT( "mie_read_word()" );
04256 return;
04257 }
04258
04259
04260 static void mie_next_data( const char *chFile,
04261 FILE *io,
04262 char chLine[],
04263 long int *dl)
04264 {
04265 char *str;
04266
04267 DEBUG_ENTRY( "mie_next_data()" );
04268
04269
04270
04271
04272
04273 chLine[0] = '#';
04274 while( chLine[0] == '#' )
04275 {
04276 mie_next_line(chFile,io,chLine,dl);
04277 }
04278
04279
04280 str = strstr(chLine,"#");
04281 if( str != NULL )
04282 *str = '\0';
04283
04284 DEBUG_EXIT( "mie_next_data()" );
04285 return;
04286 }
04287
04288
04289
04290 static void mie_next_line( const char *chFile,
04291 FILE *io,
04292 char chLine[],
04293 long int *dl)
04294 {
04295 DEBUG_ENTRY( "mie_next_line()" );
04296
04297 if( fgets( chLine, FILENAME_PATH_LENGTH_2, io ) == NULL )
04298 {
04299 fprintf( ioQQQ, " Could not read from %s\n",chFile);
04300 if( feof(io) )
04301 fprintf( ioQQQ, " EOF reached\n");
04302 puts( "[Stop in mie_next_line]" );
04303 cdEXIT(EXIT_FAILURE);
04304 }
04305 (*dl)++;
04306
04307 DEBUG_EXIT( "mie_next_line()" );
04308 return;
04309 }
04310
04311
04312
04313
04314
04315
04316
04317
04318
04319
04320
04321
04322
04323
04324 void gauss_init(long int nn,
04325 double xbot,
04326 double xtop,
04327 double x[],
04328 double a[],
04329 double rr[],
04330 double ww[])
04331 {
04332 long int i;
04333 double bma,
04334 bpa;
04335
04336 DEBUG_ENTRY( "gauss_init()" );
04337
04338 bpa = (xtop+xbot)/2.;
04339 bma = (xtop-xbot)/2.;
04340
04341 for( i=0; i < nn; i++ )
04342 {
04343 rr[i] = bpa + bma*x[nn-1-i];
04344 ww[i] = bma*a[i];
04345 }
04346
04347 DEBUG_EXIT( "gauss_init()" );
04348 return;
04349 }
04350
04351
04352
04353
04354 void gauss_legendre(long int nn,
04355 double x[],
04356 double a[])
04357 {
04358 long int i,
04359 iter,
04360 j;
04361 double *c,
04362 cc,
04363 csa,
04364 d,
04365 dp1,
04366 dpn = 0.,
04367 dq,
04368 fj,
04369 fn,
04370 pn,
04371 pn1 = 0.,
04372 q,
04373 xt = 0.;
04374
04375 const double SAFETY = 5.;
04376
04377
04378 DEBUG_ENTRY( "gauss_legendre()" );
04379
04380 if( nn%2 == 1 )
04381 {
04382 fprintf( ioQQQ, " Illegal number of abcissas\n" );
04383 puts( "[Stop in gauss_legendre]" );
04384 cdEXIT(EXIT_FAILURE);
04385 }
04386
04387 c = (double *)MALLOC(sizeof(double)*(unsigned)nn);
04388
04389 fn = (double)nn;
04390 csa = 0.;
04391 cc = 2.;
04392 for( j=1; j < nn; j++ )
04393 {
04394 fj = (double)j;
04395
04396
04397
04398 c[j] = POW2(fj)/((fj-0.5)*(fj+0.5));
04399 cc *= c[j];
04400 }
04401
04402 for( i=0; i < nn/2; i++ )
04403 {
04404 switch( i )
04405 {
04406 case 0:
04407 xt = 1. - 2.78/(4. + POW2(fn));
04408 break;
04409 case 1:
04410 xt = xt - 4.1*(1. + 0.06*(1. - 8./fn))*(1. - xt);
04411 break;
04412 case 2:
04413 xt = xt - 1.67*(1. + 0.22*(1. - 8./fn))*(x[0] - xt);
04414 break;
04415 default:
04416 xt = 3.*(x[i-1] - x[i-2]) + x[i-3];
04417 }
04418 d = 1.;
04419 for( iter=1; (iter < 20) && (fabs(d) > DBL_EPSILON); iter++ )
04420 {
04421
04422
04423
04424 pn1 = 0.5;
04425 pn = xt;
04426 dp1 = 0.;
04427 dpn = 1.;
04428 for( j=1; j < nn; j++ )
04429 {
04430
04431
04432 q = 2.*xt*pn - c[j]*pn1;
04433 dq = 2.*xt*dpn - c[j]*dp1 + 2.*pn;
04434 pn1 = pn;
04435 pn = q;
04436 dp1 = dpn;
04437 dpn = dq;
04438 }
04439 d = pn/dpn;
04440 xt -= d;
04441 }
04442 x[i] = xt;
04443 x[nn-1-i] = -xt;
04444
04445
04446 a[i] = cc/(dpn*2.*pn1);
04447 a[nn-1-i] = a[i];
04448 csa += a[i];
04449 }
04450
04451
04452
04453 if( fabs(1.-csa) > SAFETY*fn*DBL_EPSILON )
04454 {
04455 fprintf( ioQQQ, " gauss_legendre failed to converge: delta = %.4e\n",fabs(1.-csa) );
04456 puts( "[Stop in gauss_legendre]" );
04457 cdEXIT(EXIT_FAILURE);
04458 }
04459 free(c);
04460
04461 DEBUG_EXIT( "gauss_legendre()" );
04462 return;
04463 }
04464
04465
04466
04467
04468
04469
04470 void find_arr(double x,
04471 double xa[],
04472 long int n,
04473 long int *ind,
04474 bool *lgOutOfBounds)
04475 {
04476 long int i1,
04477 i2,
04478 i3,
04479 sgn,
04480 sgn2;
04481
04482 DEBUG_ENTRY( "find_arr()" );
04483
04484
04485 if( n < 2 )
04486 {
04487 fprintf( ioQQQ, " Invalid array\n");
04488 puts( "[Stop in find_arr]" );
04489 cdEXIT(EXIT_FAILURE);
04490 }
04491
04492 i1 = 0;
04493 i3 = n-1;
04494 sgn = sign3(xa[i3]-xa[i1]);
04495 if( sgn == 0 )
04496 {
04497 fprintf( ioQQQ, " Ill-ordered array\n");
04498 puts( "[Stop in find_arr]" );
04499 cdEXIT(EXIT_FAILURE);
04500 }
04501
04502 *lgOutOfBounds = x < MIN2(xa[0],xa[n-1]) || x > MAX2(xa[0],xa[n-1]);
04503 if( *lgOutOfBounds )
04504 {
04505 *ind = -1;
04506 DEBUG_EXIT( "find_arr()" );
04507 return;
04508 }
04509
04510 i2 = (n-1)/2;
04511 while( (i3-i1) > 1 )
04512 {
04513 sgn2 = sign3(x-xa[i2]);
04514 if( sgn2 != 0 )
04515 {
04516 if( sgn == sgn2 )
04517 {
04518 i1 = i2;
04519 }
04520 else
04521 {
04522 i3 = i2;
04523 }
04524 i2 = (i1+i3)/2;
04525 }
04526 else
04527 {
04528 *ind = i2;
04529 DEBUG_EXIT( "find_arr()" );
04530 return;
04531 }
04532 }
04533 *ind = i1;
04534 DEBUG_EXIT( "find_arr()" );
04535 return;
04536 }
04537
04538
04539
04540
04541
04542
04543
04544
04545
04546
04547
04548
04549
04550
04551
04552
04553
04554
04555
04556
04557
04558
04559
04560
04561
04562
04563
04564
04565 #define NMXLIM 16000
04566
04567 static void sinpar(double nre,
04568 double nim,
04569 double x,
04570 double *qext,
04571 double *qphase,
04572 double *qscat,
04573 double *ctbrqs,
04574 double *qback,
04575 long int *iflag)
04576 {
04577 long int n,
04578 nmx1,
04579 nmx2,
04580 nn,
04581 nsqbk;
04582 double ectb,
04583 eqext,
04584 eqpha,
04585 eqscat,
04586 error=0.,
04587 error1=0.,
04588 rx,
04589 t1,
04590 t2,
04591 t3,
04592 t4,
04593 t5,
04594 tx,
04595 x3,
04596 x5=0.,
04597 xcut,
04598 xrd;
04599
04600
04601
04602 complex<double> *a;
04603 complex<double> cdum1,
04604 cdum2,
04605 ci,
04606 eqb,
04607 nc,
04608 nc2,
04609 nc212,
04610 qbck,
04611 rrf,
04612 rrfx,
04613 sman,
04614 sman1,
04615 smbn,
04616 smbn1,
04617 tc1,
04618 tc2,
04619 wn,
04620 wn1,
04621 wn2;
04622
04623 DEBUG_ENTRY( "sinpar()" );
04624
04625 a = (complex<double>*)MALLOC(sizeof(complex<double>)*(unsigned)NMXLIM);
04626
04627 *iflag = 0;
04628 ci = complex<double>(0.,1.);
04629 nc = complex<double>(nre,-nim);
04630 nc2 = nc*nc;
04631 rrf = 1./nc;
04632 rx = 1./x;
04633 rrfx = rrf*rx;
04634
04635
04636
04637
04638
04639
04640
04641
04642
04643
04644
04645
04646
04647
04648
04649
04650 xrd = exp(log(x)/3.);
04651
04652
04653 t1 = x + 4.*xrd + 3.;
04654
04655
04656
04657 if( !(x <= 8. || x >= 4200.) )
04658 t1 += 0.05*xrd + 3.;
04659 t1 *= 1.01;
04660
04661
04662
04663
04664
04665
04666
04667
04668
04669 t4 = x*sqrt(nre*nre+nim*nim);
04670 nmx1 = nint(MAX2(t1,t4)) + 15;
04671
04672 if( nmx1 < NMXLIM )
04673 {
04674 nmx2 = nint(t1);
04675
04676
04677
04678
04679
04680 if( nmx2 <= 4 )
04681 {
04682 nmx2 = 4;
04683 nmx1 = nint(MAX2(4.,t4)) + 15;
04684 }
04685
04686
04687 a[nmx1] = 0.;
04688
04689
04690
04691
04692 for( n=0; n < nmx1; n++ )
04693 {
04694 nn = nmx1 - n;
04695 a[nn-1] = (double)(nn+1)*rrfx - 1./((double)(nn+1)*rrfx+a[nn]);
04696 }
04697
04698 t1 = cos(x);
04699 t2 = sin(x);
04700 wn2 = complex<double>(t1,-t2);
04701 wn1 = complex<double>(t2,t1);
04702 wn = rx*wn1 - wn2;
04703 tc1 = a[0]*rrf + rx;
04704 tc2 = a[0]*nc + rx;
04705 sman = (tc1*wn.real() - wn1.real())/(tc1*wn - wn1);
04706 smbn = (tc2*wn.real() - wn1.real())/(tc2*wn - wn1);
04707
04708
04709
04710
04711 xcut = 3.e-04;
04712 if( x < xcut )
04713 {
04714 nc212 = (nc2-1.)/(nc2+2.);
04715 x3 = POW3(x);
04716 x5 = x3*POW2(x);
04717
04718 sman = ci*2.*x3*nc212*(1./3.+x*x*0.2*(nc2-2.)/(nc2+2.)) + 4.*x5*x*nc212*nc212/9.;
04719 smbn = ci*x5*(nc2-1.)/45.;
04720 }
04721
04722 sman1 = sman;
04723 smbn1 = smbn;
04724 t1 = 1.5;
04725 sman *= t1;
04726 smbn *= t1;
04727
04728 *qext = 2.*(sman.real() + smbn.real());
04729 *qphase = 2.*(sman.imag() + smbn.imag());
04730 nsqbk = -1;
04731 qbck = -2.*(sman - smbn);
04732 *qscat = (norm(sman) + norm(smbn))/.75;
04733
04734 *ctbrqs = 0.0;
04735 n = 2;
04736
04737
04738 while( true )
04739 {
04740 t1 = 2.*(double)n - 1.;
04741 t3 = 2.*(double)n + 1.;
04742 wn2 = wn1;
04743 wn1 = wn;
04744 wn = t1*rx*wn1 - wn2;
04745 cdum1 = a[n-1];
04746 cdum2 = n*rx;
04747 tc1 = cdum1*rrf + cdum2;
04748 tc2 = cdum1*nc + cdum2;
04749 sman = (tc1*wn.real() - wn1.real())/(tc1*wn - wn1);
04750 smbn = (tc2*wn.real() - wn1.real())/(tc2*wn - wn1);
04751
04752
04753
04754 if( x < xcut && n == 2 )
04755 {
04756
04757 sman = ci*x5*(nc2-1.)/(15.*(2.*nc2+3.));
04758 smbn = 0.;
04759 }
04760
04761 eqext = t3*(sman.real() + smbn.real());
04762 *qext += eqext;
04763 eqpha = t3*(sman.imag() + smbn.imag());
04764 *qphase += eqpha;
04765 nsqbk = -nsqbk;
04766 eqb = t3*(sman - smbn)*(double)nsqbk;
04767 qbck += eqb;
04768 tx = norm(sman) + norm(smbn);
04769 eqscat = t3*tx;
04770 *qscat += eqscat;
04771 t2 = (double)(n - 1);
04772 t5 = (double)n;
04773 t4 = t1/(t5*t2);
04774 t2 = (t2*(t5 + 1.))/t5;
04775 ectb = t2*(sman1.real()*sman.real()+sman1.imag()*sman.imag() + smbn1.real()*smbn.real() +
04776 smbn1.imag()*smbn.imag()) +
04777 t4*(sman1.real()*smbn1.real()+sman1.imag()*smbn1.imag());
04778 *ctbrqs += ectb;
04779
04780
04781
04782 if( tx < 1.e-14 )
04783 {
04784
04785 eqext = fabs(eqext/ *qext);
04786 eqpha = fabs(eqpha/ *qphase);
04787 eqscat = fabs(eqscat/ *qscat);
04788 ectb = ( n == 2 ) ? 0. : fabs(ectb/ *ctbrqs);
04789 eqb = complex<double>( fabs(eqb.real()/qbck.real()), fabs(eqb.imag()/qbck.imag()) );
04790
04791 error = MAX4(eqext,eqpha,eqscat,ectb);
04792
04793 error1 = MAX2(eqb.real(),eqb.imag());
04794 if( error < 1.e-07 && error1 < 1.e-04 )
04795 break;
04796
04797
04798
04799
04800 if( x < xcut )
04801 break;
04802 }
04803
04804 smbn1 = smbn;
04805 sman1 = sman;
04806 n++;
04807 if( n > nmx2 )
04808 {
04809 *iflag = 1;
04810 break;
04811 }
04812 }
04813
04814 t1 = 2.*POW2(rx);
04815 *qext *= t1;
04816 *qphase *= t1;
04817 *qback = norm(qbck)*POW2(rx);
04818 *qscat *= t1;
04819 *ctbrqs *= 2.*t1;
04820 }
04821 else
04822 {
04823 *iflag = 2;
04824 }
04825
04826 free( a );
04827
04828 DEBUG_EXIT( "sinpar()" );
04829 return;
04830 }
04831
04832
04833 static void anomal(double x,
04834 double *qext,
04835 double *qabs,
04836 double *qphase,
04837 double *xistar,
04838 double delta,
04839 double beta)
04840 {
04841
04842
04843
04844
04845
04846
04847
04848 double xi,
04849 xii;
04850 complex<double> cbigk,
04851 ci,
04852 cw;
04853
04854 DEBUG_ENTRY( "anomal()" );
04855
04856
04857 xi = 2.*x*delta;
04858 xii = 2.*x*beta;
04859
04860 *xistar = sqrt(POW2(xi)+POW2(xii));
04861
04862 ci = complex<double>(0.,1.);
04863 cw = -complex<double>(xi,xii)*ci;
04864 bigk(cw,&cbigk);
04865 *qext = 4.*cbigk.real();
04866 *qphase = 4.*cbigk.imag();
04867 cw = 2.*xii;
04868 bigk(cw,&cbigk);
04869 *qabs = 2.*cbigk.real();
04870
04871
04872 DEBUG_EXIT( "anomal()" );
04873 return;
04874 }
04875
04876
04877 static void bigk(complex<double> cw,
04878 complex<double> *cbigk)
04879 {
04880
04881
04882
04883
04884 DEBUG_ENTRY( "bigk()" );
04885
04886 if( abs(cw) < 1.e-2 )
04887 {
04888
04889
04890
04891 *cbigk = cw*((1./3.)-cw*((1./8.)-cw*((1./30.)-cw*((1./144.)-cw*((1./840.)-cw*(1./5760.))))));
04892 }
04893 else
04894 {
04895 *cbigk = 0.5 + (exp(-cw)*(1.+cw)-1.)/(cw*cw);
04896 }
04897
04898 DEBUG_EXIT( "bigk()" );
04899 return;
04900 }
04901
04902
04903
04904 static void ritodf(double nr,
04905 double ni,
04906 double *eps1,
04907 double *eps2)
04908 {
04909
04910 DEBUG_ENTRY( "ritodf()" );
04911
04912 *eps1 = nr*nr - ni*ni;
04913 *eps2 = 2.*nr*ni;
04914
04915 DEBUG_EXIT( "ritodf()" );
04916 return;
04917 }
04918
04919
04920
04921 static void dftori( double *nr,
04922 double *ni,
04923 double eps1,
04924 double eps2)
04925 {
04926 double eps;
04927
04928 DEBUG_ENTRY( "dftori()" );
04929
04930 eps = sqrt(eps2*eps2+eps1*eps1);
04931 *nr = sqrt((eps+eps1)/2.);
04932 ASSERT( *nr > 0. );
04933
04934
04935
04936 *ni = eps2/(2.*(*nr));
04937
04938 DEBUG_EXIT( "dftori()" );
04939 return;
04940 }
04941