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