00001
00002
00003
00004
00005
00006
00007
00008
00009 #include "cddefines.h"
00010 #include "rfield.h"
00011 #include "iterations.h"
00012 #include "physconst.h"
00013 #include "doppvel.h"
00014 #include "dense.h"
00015 #include "trace.h"
00016 #include "opacity.h"
00017 #include "ipoint.h"
00018 #include "geometry.h"
00019 #include "continuum.h"
00020
00021
00022 STATIC void read_continuum_mesh( void );
00023
00024
00025 STATIC void fill(double fenlo,
00026 double fenhi,
00027 double resolv,
00028 long int *n0,
00029 long int *ipnt,
00030
00031 bool lgCount );
00032
00033
00034 STATIC void rfield_opac_malloc(void);
00035
00036
00037 STATIC void ChckFill(void);
00038
00039 void ContCreateMesh(void)
00040 {
00041 long int
00042 i,
00043 ipnt,
00044 n0;
00045
00046
00047 static bool lgPntEval = false;
00048
00049 DEBUG_ENTRY( "ContCreateMesh()" );
00050
00051
00052
00053
00054 if( lgPntEval )
00055 {
00056 if( trace.lgTrace )
00057 {
00058 fprintf( ioQQQ, " ContCreateMesh called, not evaluating.\n" );
00059 }
00060
00061 for( i=0; i < rfield.nupper; i++ )
00062 {
00063 rfield.anu[i] = rfield.AnuOrg[i];
00064 rfield.anu2[i] = rfield.anu[i]*rfield.anu[i];
00065 }
00066 return;
00067 }
00068 else
00069 {
00070 if( trace.lgTrace )
00071 {
00072 fprintf( ioQQQ, " ContCreateMesh called first time.\n" );
00073 }
00074 lgPntEval = true;
00075 }
00076
00077
00078
00079
00080
00081 read_continuum_mesh();
00082
00083
00084
00085
00086
00087
00088
00089
00090
00091
00092
00093
00094
00095
00096 n0 = 1;
00097 ipnt = 0;
00098
00099 continuum.nrange = 0;
00100
00101
00102
00103 n0 = 1;
00104 ipnt = 0;
00105
00106 continuum.nrange = 0;
00107 continuum.StoredEnergy[continuum.nStoredBands-1] = rfield.egamry;
00108
00109 fill(rfield.emm, continuum.StoredEnergy[0] , continuum.StoredResolution[0],&n0,&ipnt,true );
00110 for(i=1; i<continuum.nStoredBands; ++i )
00111 {
00112 fill(continuum.StoredEnergy[i-1] ,
00113 continuum.StoredEnergy[i],
00114 continuum.StoredResolution[i],
00115 &n0,&ipnt,true);
00116 }
00117
00118
00119
00120
00121 rfield.nupper = n0 - 1;
00122
00123 if( rfield.nupper >= NCELL )
00124 {
00125 fprintf(ioQQQ," Currently the arrays that hold interpolated tables can only hold %i points.\n",NCELL);
00126 fprintf(ioQQQ," This continuum mesh really needs to have %li points.\n",rfield.nupper);
00127 fprintf(ioQQQ," Please increase the value of NCELL in rfield.h and recompile.\n Sorry.");
00128 cdEXIT(EXIT_FAILURE);
00129 }
00130
00131
00132
00133 rfield.nflux = rfield.nupper-1;
00134
00135
00136
00137 rfield_opac_malloc();
00138
00139
00140
00141
00142 geometry.nend_max = geometry.nend[0];
00143 for( i=1; i < iterations.iter_malloc; i++ )
00144 {
00145 geometry.nend_max = MAX2( geometry.nend_max , geometry.nend[i] );
00146 }
00147
00148 rfield.ConEmitLocal = (realnum**)MALLOC( (size_t)(geometry.nend_max+1)*sizeof(realnum *) );
00149 for( i=0;i<(geometry.nend_max+1); ++i )
00150 {
00151 rfield.ConEmitLocal[i] = (realnum*)MALLOC( (size_t)rfield.nupper*sizeof(realnum) );
00152 }
00153
00154
00155 n0 = 1;
00156 ipnt = 0;
00157
00158
00159 continuum.nrange = 0;
00160
00161
00162 fill(rfield.emm, continuum.StoredEnergy[0] , continuum.StoredResolution[0] ,
00163 &n0,&ipnt,false);
00164 for(i=1; i<continuum.nStoredBands; ++i )
00165 {
00166 fill(continuum.StoredEnergy[i-1] ,
00167 continuum.StoredEnergy[i],
00168 continuum.StoredResolution[i],
00169 &n0,&ipnt,false);
00170 }
00171
00172
00173
00174
00175 rfield.widflx[rfield.nupper] = rfield.widflx[rfield.nupper-1];
00176 rfield.anu[rfield.nupper] = rfield.anu[rfield.nupper-1] +
00177 rfield.widflx[rfield.nupper];
00178
00179
00180
00181
00182
00183
00184
00185 rfield_opac_zero( 0 , rfield.nupper );
00186
00187
00188 ChckFill();
00189
00190
00191 for( i=1; i<rfield.nupper-1; ++i )
00192 {
00193 rfield.widflx[i] = ((rfield.anu[i+1] - rfield.anu[i]) + (rfield.anu[i] -
00194 rfield.anu[i-1]))/2.f;
00195 }
00196
00197 ipnt = 0;
00198
00199 for( i=0; i < rfield.nupper; i++ )
00200 {
00201 double alf , bet;
00202
00203 rfield.AnuOrg[i] = rfield.anu[i];
00204 rfield.anusqr[i] = (realnum)sqrt(rfield.AnuOrg[i]);
00205
00206
00207 alf = 1./(1. + rfield.anu[i]*(1.1792e-4 + 7.084e-10*rfield.anu[i]));
00208 bet = 1. - alf*rfield.anu[i]*(1.1792e-4 + 2.*7.084e-10*rfield.anu[i])/4.;
00209 rfield.csigh[i] = (realnum)(alf*rfield.anu[i]*rfield.anu[i]*3.858e-25);
00210 rfield.csigc[i] = (realnum)(alf*bet*rfield.anu[i]*3.858e-25);
00211 rfield.anu2[i] = rfield.anu[i]*rfield.anu[i];
00212
00213
00214
00215 if( rfield.anu[i] < rfield.fine_ener_lo || rfield.anu[i]>rfield.fine_ener_hi )
00216 {
00217
00218 rfield.ipnt_coarse_2_fine[i] = 0;
00219 }
00220 else
00221 {
00222 if( ipnt==0 )
00223 {
00224
00225 rfield.ipnt_coarse_2_fine[i] = 0;
00226 ipnt = 1;
00227 }
00228 else
00229 {
00230
00231 while (ipnt < rfield.nfine_malloc && rfield.fine_anu[ipnt]<rfield.anu[i] )
00232 {
00233 ++ipnt;
00234 }
00235 rfield.ipnt_coarse_2_fine[i] = ipnt;
00236 }
00237 }
00238
00239
00240 rfield.trans_coef_zone[i] = 1.f;
00241 rfield.trans_coef_total[i] = 1.f;
00242 }
00243 return;
00244 }
00245
00246
00247 STATIC void fill(
00248
00249 double fenlo,
00250
00251 double fenhi,
00252
00253 double resolv,
00254
00255 long int *n0,
00256
00257 long int *ipnt,
00258
00259 bool lgCount )
00260 {
00261 long int i,
00262 nbin;
00263 realnum widtot;
00264 double aaa , bbb;
00265
00266 DEBUG_ENTRY( "fill()" );
00267
00268 ASSERT( fenlo>0. && fenhi>0. && resolv>0. );
00269
00270
00271 nbin = (long int)(log(10.)*log10(fenhi/fenlo)/resolv + 1);
00272
00273 if( lgCount )
00274 {
00275
00276 *n0 += nbin;
00277 return;
00278 }
00279
00280 if( *ipnt > 0 && fabs(1.-fenlo/continuum.filbnd[*ipnt]) > 1e-4 )
00281 {
00282 fprintf( ioQQQ, " FILL improper bounds.\n" );
00283 fprintf( ioQQQ, " ipnt=%3ld fenlo=%11.4e filbnd(ipnt)=%11.4e\n",
00284 *ipnt, fenlo, continuum.filbnd[*ipnt] );
00285 cdEXIT(EXIT_FAILURE);
00286 }
00287
00288 ASSERT( *ipnt < continuum.nStoredBands );
00289
00290 continuum.ifill0[*ipnt] = *n0 - 1;
00291 continuum.filbnd[*ipnt] = (realnum)fenlo;
00292 continuum.filbnd[*ipnt+1] = (realnum)fenhi;
00293
00294
00295
00296 continuum.fildel[*ipnt] = (realnum)(log10(fenhi/fenlo)/nbin);
00297
00298 if( continuum.fildel[*ipnt] < 0.01 )
00299 {
00300 continuum.filres[*ipnt] = (realnum)(log(10.)*continuum.fildel[*ipnt]);
00301 }
00302 else
00303 {
00304 continuum.filres[*ipnt] = (realnum)((pow(10.,2.*continuum.fildel[*ipnt]) - 1.)/2./
00305 pow((realnum)10.f,continuum.fildel[*ipnt]));
00306 }
00307
00308 if( (*n0 + nbin-2) > rfield.nupper )
00309 {
00310 fprintf( ioQQQ, " Fill would need %ld cells to get to an energy of %.3e\n",
00311 *n0 + nbin, fenhi );
00312 fprintf( ioQQQ, " This is a major logical error in fill.\n");
00313 ShowMe();
00314 cdEXIT(EXIT_FAILURE);
00315 }
00316
00317 widtot = 0.;
00318 for( i=0; i < nbin; i++ )
00319 {
00320 bbb = continuum.fildel[*ipnt]*((realnum)(i) + 0.5);
00321 aaa = pow( 10. , bbb );
00322
00323 rfield.anu[i+continuum.ifill0[*ipnt]] = (realnum)(fenlo*aaa);
00324
00325 rfield.widflx[i+continuum.ifill0[*ipnt]] = rfield.anu[i+continuum.ifill0[*ipnt]]*
00326 continuum.filres[*ipnt];
00327
00328 widtot += rfield.widflx[i+continuum.ifill0[*ipnt]];
00329 }
00330
00331 *n0 += nbin;
00332 if( trace.lgTrace && (trace.lgConBug || trace.lgPtrace) )
00333 {
00334 fprintf( ioQQQ,
00335 " FILL range%2ld from%10.3e to%10.3eR in%4ld cell; ener res=%10.3e WIDTOT=%10.3e\n",
00336 *ipnt,
00337 rfield.anu[continuum.ifill0[*ipnt]] - rfield.widflx[continuum.ifill0[*ipnt]]/2.,
00338 rfield.anu[continuum.ifill0[*ipnt]+nbin-1] + rfield.widflx[continuum.ifill0[*ipnt]+nbin-1]/2.,
00339 nbin,
00340 continuum.filres[*ipnt],
00341 widtot );
00342
00343 fprintf( ioQQQ, " The requested range was%10.3e%10.3e The requested resolution was%10.3e\n",
00344 fenlo, fenhi, resolv );
00345 }
00346
00347
00348 *ipnt += 1;
00349 continuum.nrange = MAX2(continuum.nrange,*ipnt);
00350 return;
00351 }
00352
00353
00354 STATIC void ChckFill(void)
00355 {
00356 bool lgFail;
00357 long int i,
00358 ipnt;
00359 double energy;
00360
00361 DEBUG_ENTRY( "ChckFill()" );
00362
00363 ASSERT( rfield.anu[0] >= rfield.emm*0.99 );
00364 ASSERT( rfield.anu[rfield.nupper-1] <= rfield.egamry*1.01 );
00365
00366 lgFail = false;
00367 for( i=0; i < continuum.nrange; i++ )
00368 {
00369
00370 energy = (continuum.filbnd[i] + continuum.filbnd[i+1])/2.;
00371 ipnt = ipoint(energy);
00372 if( energy < rfield.anu[ipnt-1] - rfield.widflx[ipnt-1]*0.5 )
00373 {
00374 fprintf( ioQQQ, " ChckFill middle test low fail\n" );
00375 lgFail = true;
00376 }
00377
00378
00379
00380
00381 else if( (energy > rfield.anu[ipnt-1] + rfield.widflx[ipnt-1]*0.5) &&
00382 ( energy > rfield.anu[ipnt] - rfield.widflx[ipnt]*0.5 ) )
00383 {
00384 fprintf( ioQQQ, " ChckFill middle test high fail\n" );
00385 lgFail = true;
00386 }
00387
00388
00389 energy = continuum.filbnd[i]*0.99 + continuum.filbnd[i+1]*0.01;
00390 ipnt = ipoint(energy);
00391 if( energy < rfield.anu[ipnt-1] - rfield.widflx[ipnt-1]*0.5 )
00392 {
00393 fprintf( ioQQQ, " ChckFill low test low fail\n" );
00394 lgFail = true;
00395 }
00396
00397 else if( energy > rfield.anu[ipnt-1] + rfield.widflx[ipnt-1]* 0.5 )
00398 {
00399 fprintf( ioQQQ, " ChckFill low test high fail\n" );
00400 lgFail = true;
00401 }
00402
00403
00404 energy = continuum.filbnd[i]*0.01 + continuum.filbnd[i+1]*0.99;
00405 ipnt = ipoint(energy);
00406
00407 if( energy < rfield.anu[ipnt-1] - rfield.widflx[ipnt-1]*0.5 )
00408 {
00409 fprintf( ioQQQ, " ChckFill high test low fail\n" );
00410 lgFail = true;
00411 }
00412
00413
00414
00415 else if( (energy > rfield.anu[ipnt-1] + rfield.widflx[ipnt-1]*0.5) &&
00416 ( energy > rfield.anu[ipnt] - rfield.widflx[ipnt]*0.5 ) )
00417 {
00418 fprintf( ioQQQ, " ChckFill high test high fail\n" );
00419 lgFail = true;
00420 }
00421 }
00422
00423 if( lgFail )
00424 {
00425 cdEXIT(EXIT_FAILURE);
00426 }
00427 return;
00428 }
00429
00430
00431 STATIC void rfield_opac_malloc(void)
00432 {
00433 long i;
00434
00435 DEBUG_ENTRY( "rfield_opac_malloc()" );
00436
00437
00438
00439 ++rfield.nupper;
00440
00441
00446
00447 rfield.fine_ener_lo = 0.05f;
00448 rfield.fine_ener_hi = 1500.f;
00449
00450
00451
00452
00453
00454
00455
00456
00457
00458
00459
00460
00461
00462
00463
00466 double TeLowestFineOpacity = 1e4;
00467 rfield.fine_opac_velocity_width =
00468 (realnum)sqrt(2.*BOLTZMANN/ATOMIC_MASS_UNIT*TeLowestFineOpacity/
00469 dense.AtomicWeight[rfield.fine_opac_nelem] ) /
00470
00471
00472 rfield.fine_opac_nresolv;
00473
00474
00475 rfield.ipFineConVelShift = 0;
00476
00477
00478 rfield.fine_resol = rfield.fine_opac_velocity_width / SPEEDLIGHT;
00479
00480
00481 rfield.nfine_malloc = (long)(log10( rfield.fine_ener_hi / rfield.fine_ener_lo ) / log10( 1. + rfield.fine_resol ) );
00482 if( rfield.nfine_malloc <= 0 )
00483 TotalInsanity();
00484 rfield.nfine = rfield.nfine_malloc;
00485
00486
00487 rfield.fine_opac_zone = (realnum *)MALLOC(sizeof(realnum)*(unsigned)rfield.nfine_malloc );
00488 memset(rfield.fine_opac_zone , 0 , (unsigned long)rfield.nfine_malloc*sizeof(realnum) );
00489
00490
00491 rfield.fine_opt_depth = (realnum *)MALLOC(sizeof(realnum)*(unsigned)rfield.nfine_malloc );
00492 memset(rfield.fine_opt_depth , 0 , (unsigned long)rfield.nfine_malloc*sizeof(realnum) );
00493
00494 rfield.fine_anu = (realnum *)MALLOC(sizeof(realnum)*(unsigned)rfield.nfine_malloc );
00495
00496
00497 ASSERT( rfield.fine_ener_lo > 0. && rfield.fine_resol > 0 );
00498 for( i=0;i<rfield.nfine_malloc; ++i )
00499 {
00500 rfield.fine_anu[i] = rfield.fine_ener_lo * (realnum)pow( (1.+rfield.fine_resol), (i+1.) );
00501 }
00502
00503
00504
00505 rfield.line_count = (long *)MALLOC(sizeof(long)*(unsigned)NCELL );
00506 for( i=0; i<rfield.nupper; ++i)
00507 {
00508 rfield.line_count[i] = 0;
00509 }
00510 rfield.anu = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
00511 rfield.AnuOrg = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
00512 rfield.widflx = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
00513 rfield.anulog = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
00514 rfield.anusqr = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
00515 rfield.anu2 = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
00516 rfield.anu3 = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
00517 rfield.flux_beam_time = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
00518 rfield.flux_isotropic = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
00519 rfield.flux_beam_const = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
00520 rfield.flux = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
00521 rfield.flux_accum = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
00522 rfield.convoc = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
00523 rfield.OccNumbBremsCont = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
00524 rfield.OccNumbIncidCont = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
00525 rfield.OccNumbDiffCont = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
00526 rfield.OccNumbContEmitOut = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
00527 rfield.ConEmitReflec = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
00528 rfield.ConEmitOut = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
00529 rfield.ConInterOut = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
00530 rfield.ConRefIncid = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
00531 rfield.SummedCon = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
00532 rfield.SummedDif = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
00533 rfield.SummedDifSave = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
00534 rfield.SummedOcc = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
00535 rfield.ConOTS_local_photons = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
00536 rfield.TotDiff2Pht = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
00537 rfield.ConOTS_local_OTS_rate = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
00538 rfield.otslin = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
00539 rfield.otscon = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
00540 rfield.otslinNew = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
00541 rfield.otsconNew = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
00542 rfield.outlin = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
00543 rfield.outlin_noplot = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
00544 rfield.reflin = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
00545 rfield.flux_total_incident = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
00546 rfield.flux_beam_const_save = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
00547 rfield.flux_time_beam_save = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
00548 rfield.flux_isotropic_save = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
00549 rfield.trans_coef_zone = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
00550 rfield.trans_coef_total = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
00551 rfield.ipnt_coarse_2_fine = (long int*)MALLOC((size_t)(rfield.nupper*sizeof(long int)) );
00552
00553
00554
00555 rfield.gff = (realnum**)MALLOC((size_t)((LIMELM+1)*sizeof(realnum*)) );
00556 for( i = 1; i <= LIMELM; i++ )
00557 {
00558 rfield.gff[i] = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
00559 }
00560
00561 rfield.csigh = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
00562 rfield.csigc = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
00563 rfield.ConTabRead = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
00564
00565 rfield.comdn = (double*)MALLOC((size_t)(rfield.nupper*sizeof(double)) );
00566 rfield.comup = (double*)MALLOC((size_t)(rfield.nupper*sizeof(double)) );
00567 rfield.ContBoltz = (double*)MALLOC((size_t)(rfield.nupper*sizeof(double)) );
00568
00569
00570 rfield.otssav = (realnum**)MALLOC((size_t)(rfield.nupper*sizeof(realnum*)));
00571 for( i=0; i<rfield.nupper; ++i)
00572 {
00573 rfield.otssav[i] = (realnum*)MALLOC(2*sizeof(realnum));
00574 }
00575
00576
00577
00578 rfield.chLineLabel = (char**)MALLOC((size_t)(rfield.nupper*sizeof(char*)));
00579 rfield.chContLabel = (char**)MALLOC((size_t)(rfield.nupper*sizeof(char*)));
00580
00581
00582 for( i=0; i<rfield.nupper; ++i)
00583 {
00584 rfield.chLineLabel[i] = (char*)MALLOC(5*sizeof(char));
00585 rfield.chContLabel[i] = (char*)MALLOC(5*sizeof(char));
00586 }
00587
00588 opac.TauAbsFace = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
00589 memset( opac.TauAbsFace , 0 , rfield.nupper*sizeof(realnum) );
00590
00591 opac.TauScatFace = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
00592 opac.E2TauAbsFace = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
00593 opac.E2TauAbsTotal = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
00594 opac.TauAbsTotal = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
00595 opac.E2TauAbsOut = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
00596 opac.ExpmTau = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
00597 opac.tmn = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
00598
00599 opac.opacity_abs = (double*)MALLOC((size_t)(rfield.nupper*sizeof(double)) );
00600 opac.opacity_abs_savzon1 = (double*)MALLOC((size_t)(rfield.nupper*sizeof(double)) );
00601 opac.OldOpacSave = (double*)MALLOC((size_t)(rfield.nupper*sizeof(double)) );
00602 opac.opacity_sct = (double*)MALLOC((size_t)(rfield.nupper*sizeof(double)) );
00603 opac.albedo = (double*)MALLOC((size_t)(rfield.nupper*sizeof(double)) );
00604 opac.opacity_sct_savzon1 = (double*)MALLOC((size_t)(rfield.nupper*sizeof(double)) );
00605 opac.OpacStatic = (double*)MALLOC((size_t)(rfield.nupper*sizeof(double)) );
00606 opac.FreeFreeOpacity = (double*)MALLOC((size_t)(rfield.nupper*sizeof(double)) );
00607 opac.ExpZone = (double*)MALLOC((size_t)(rfield.nupper*sizeof(double)) );
00608
00609 opac.TauAbsGeo = (realnum**)MALLOC((size_t)(2*sizeof(realnum *)) );
00610 opac.TauScatGeo = (realnum**)MALLOC((size_t)(2*sizeof(realnum *)) );
00611 opac.TauTotalGeo = (realnum**)MALLOC((size_t)(2*sizeof(realnum *)) );
00612
00613 for( i=0; i<2; ++i)
00614 {
00615 opac.TauAbsGeo[i] = (realnum*)MALLOC(rfield.nupper*sizeof(realnum));
00616 opac.TauScatGeo[i] = (realnum*)MALLOC(rfield.nupper*sizeof(realnum));
00617 opac.TauTotalGeo[i] = (realnum*)MALLOC(rfield.nupper*sizeof(realnum));
00618 }
00619
00620
00621 --rfield.nupper;
00622
00623
00624 lgRfieldMalloced = true;
00625 return;
00626 }
00627
00628
00629
00630 STATIC void read_continuum_mesh( void )
00631 {
00632 FILE *ioDATA;
00633 char chLine[INPUT_LINE_LENGTH];
00634 long i;
00635 bool lgEOL;
00636 long i1 , i2 , i3;
00637
00638 DEBUG_ENTRY( "read_continuum_mesh()" );
00639
00640 if( trace.lgTrace )
00641 fprintf( ioQQQ," read_continuum_mesh opening continuum_mesh.ini:");
00642
00643 ioDATA = open_data( "continuum_mesh.ini", "r" );
00644
00645
00646 if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
00647 {
00648 fprintf( ioQQQ, " read_continuum_mesh could not read first line of continuum_mesh.ini.\n");
00649 cdEXIT(EXIT_FAILURE);
00650 }
00651
00652
00653 continuum.nStoredBands = 0;
00654 while( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) != NULL )
00655 {
00656
00657
00658 if( chLine[0] != '#')
00659 ++continuum.nStoredBands;
00660 }
00661
00662
00663
00664 continuum.filbnd =
00665 ((realnum *)MALLOC( (size_t)(continuum.nStoredBands+1)*sizeof(realnum )));
00666 continuum.fildel =
00667 ((realnum *)MALLOC( (size_t)(continuum.nStoredBands+1)*sizeof(realnum )));
00668 continuum.filres =
00669 ((realnum *)MALLOC( (size_t)(continuum.nStoredBands+1)*sizeof(realnum )));
00670 continuum.ifill0 =
00671 ((long *)MALLOC( (size_t)(continuum.nStoredBands+1)*sizeof(long )));
00672 continuum.StoredEnergy =
00673 ((double *)MALLOC( (size_t)(continuum.nStoredBands+1)*sizeof(double )));
00674 continuum.StoredResolution =
00675 ((double *)MALLOC( (size_t)(continuum.nStoredBands+1)*sizeof(double )));
00676
00677
00678 if( fseek( ioDATA , 0 , SEEK_SET ) != 0 )
00679 {
00680 fprintf( ioQQQ, " read_continuum_mesh could not rewind continuum_mesh.ini.\n");
00681 cdEXIT(EXIT_FAILURE);
00682 }
00683
00684
00685 if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
00686 {
00687 fprintf( ioQQQ, " read_continuum_mesh could not read first line of continuum_mesh.ini.\n");
00688 cdEXIT(EXIT_FAILURE);
00689 }
00690
00691 i = 1;
00692
00693 i1 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
00694 i2 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
00695 i3 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
00696
00697
00698 if( ( i1 != 1 ) || ( i2 != 9 ) || ( i3 != 29 ) )
00699 {
00700 fprintf( ioQQQ,
00701 " read_continuum_mesh: the version of continuum_mesh.ini is not the current version.\n" );
00702 fprintf( ioQQQ,
00703 " I expected to find the number 01 09 29 and got %li %li %li instead.\n" ,
00704 i1 , i2 , i3 );
00705 fprintf( ioQQQ, "Here is the line image:\n==%s==\n", chLine );
00706 cdEXIT(EXIT_FAILURE);
00707 }
00708
00709
00710
00711 continuum.nStoredBands = 0;
00712 while( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) != NULL )
00713 {
00714
00715 if( chLine[0] != '#')
00716 {
00717 i = 1;
00718 continuum.StoredEnergy[continuum.nStoredBands] = FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
00719 continuum.StoredResolution[continuum.nStoredBands] = FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
00720
00721
00722 if( continuum.StoredEnergy[continuum.nStoredBands]<0. )
00723 continuum.StoredEnergy[continuum.nStoredBands] =
00724 pow(10.,continuum.StoredEnergy[continuum.nStoredBands]);
00725
00726 if( continuum.StoredResolution[continuum.nStoredBands]<0. )
00727 continuum.StoredResolution[continuum.nStoredBands] =
00728 pow(10.,continuum.StoredResolution[continuum.nStoredBands]);
00729
00730
00731 continuum.StoredResolution[continuum.nStoredBands] *= continuum.ResolutionScaleFactor;
00732
00733 if( continuum.StoredResolution[continuum.nStoredBands] == 0. )
00734 {
00735 fprintf( ioQQQ,
00736 " read_continuum_mesh: A continuum resolution was zero - this is not allowed.\n" );
00737 cdEXIT(EXIT_FAILURE);
00738 }
00739
00740 ++continuum.nStoredBands;
00741 }
00742 }
00743
00744 fclose( ioDATA );
00745
00746
00747 for( i=1; i<continuum.nStoredBands-1; ++i )
00748 {
00749 if( continuum.StoredEnergy[i-1]>=continuum.StoredEnergy[i] )
00750 {
00751 fprintf( ioQQQ,
00752 " read_continuum_mesh: The continuum definition array energies must be in increasing order.\n" );
00753 cdEXIT(EXIT_FAILURE);
00754 }
00755 }
00756 if( continuum.StoredEnergy[continuum.nStoredBands-1]!=0 )
00757 {
00758 fprintf( ioQQQ,
00759 " read_continuum_mesh: The last continuum array energies must be zero.\n" );
00760 cdEXIT(EXIT_FAILURE);
00761 }
00762 return;
00763 }