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