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 rfield.trans_coef_zone[i] = 1.f;
00249 rfield.trans_coef_total[i] = 1.f;
00250 }
00251 return;
00252 }
00253
00254
00255 STATIC void fill(
00256
00257 double fenlo,
00258
00259 double fenhi,
00260
00261 double resolv,
00262
00263 long int *n0,
00264
00265 long int *ipnt,
00266
00267 bool lgCount )
00268 {
00269 long int i,
00270 nbin;
00271 realnum widtot;
00272 double aaa , bbb;
00273
00274 DEBUG_ENTRY( "fill()" );
00275
00276 ASSERT( fenlo>0. && fenhi>0. && resolv>0. );
00277
00278
00279 nbin = (long int)(log(10.)*log10(fenhi/fenlo)/resolv + 1);
00280
00281 if( lgCount )
00282 {
00283
00284 *n0 += nbin;
00285 return;
00286 }
00287
00288 if( *ipnt > 0 && fabs(1.-fenlo/continuum.filbnd[*ipnt]) > 1e-4 )
00289 {
00290 fprintf( ioQQQ, " FILL improper bounds.\n" );
00291 fprintf( ioQQQ, " ipnt=%3ld fenlo=%11.4e filbnd(ipnt)=%11.4e\n",
00292 *ipnt, fenlo, continuum.filbnd[*ipnt] );
00293 cdEXIT(EXIT_FAILURE);
00294 }
00295
00296 ASSERT( *ipnt < continuum.nStoredBands );
00297
00298 continuum.ifill0[*ipnt] = *n0 - 1;
00299 continuum.filbnd[*ipnt] = (realnum)fenlo;
00300 continuum.filbnd[*ipnt+1] = (realnum)fenhi;
00301
00302
00303
00304 continuum.fildel[*ipnt] = (realnum)(log10(fenhi/fenlo)/nbin);
00305
00306 if( continuum.fildel[*ipnt] < 0.01 )
00307 {
00308 continuum.filres[*ipnt] = (realnum)(log(10.)*continuum.fildel[*ipnt]);
00309 }
00310 else
00311 {
00312 continuum.filres[*ipnt] = (realnum)((pow(10.,2.*continuum.fildel[*ipnt]) - 1.)/2./
00313 pow((realnum)10.f,continuum.fildel[*ipnt]));
00314 }
00315
00316 if( (*n0 + nbin-2) > rfield.nupper )
00317 {
00318 fprintf( ioQQQ, " Fill would need %ld cells to get to an energy of %.3e\n",
00319 *n0 + nbin, fenhi );
00320 fprintf( ioQQQ, " This is a major logical error in fill.\n");
00321 ShowMe();
00322 cdEXIT(EXIT_FAILURE);
00323 }
00324
00325 widtot = 0.;
00326 for( i=0; i < nbin; i++ )
00327 {
00328 bbb = continuum.fildel[*ipnt]*((realnum)(i) + 0.5);
00329 aaa = pow( 10. , bbb );
00330
00331 rfield.anu[i+continuum.ifill0[*ipnt]] = (realnum)(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 return;
00359 }
00360
00361
00362 STATIC void ChckFill(void)
00363 {
00364 bool lgFail;
00365 long int i,
00366 ipnt;
00367 double energy;
00368
00369 DEBUG_ENTRY( "ChckFill()" );
00370
00371 ASSERT( rfield.anu[0] >= rfield.emm*0.99 );
00372 ASSERT( rfield.anu[rfield.nupper-1] <= rfield.egamry*1.01 );
00373
00374 lgFail = false;
00375 for( i=0; i < continuum.nrange; i++ )
00376 {
00377
00378 energy = (continuum.filbnd[i] + continuum.filbnd[i+1])/2.;
00379 ipnt = ipoint(energy);
00380 if( energy < rfield.anu[ipnt-1] - rfield.widflx[ipnt-1]*0.5 )
00381 {
00382 fprintf( ioQQQ, " ChckFill middle test low fail\n" );
00383 lgFail = true;
00384 }
00385
00386
00387
00388
00389 else if( (energy > rfield.anu[ipnt-1] + rfield.widflx[ipnt-1]*0.5) &&
00390 ( energy > rfield.anu[ipnt] - rfield.widflx[ipnt]*0.5 ) )
00391 {
00392 fprintf( ioQQQ, " ChckFill middle test high fail\n" );
00393 lgFail = true;
00394 }
00395
00396
00397 energy = continuum.filbnd[i]*0.99 + continuum.filbnd[i+1]*0.01;
00398 ipnt = ipoint(energy);
00399 if( energy < rfield.anu[ipnt-1] - rfield.widflx[ipnt-1]*0.5 )
00400 {
00401 fprintf( ioQQQ, " ChckFill low test low fail\n" );
00402 lgFail = true;
00403 }
00404
00405 else if( energy > rfield.anu[ipnt-1] + rfield.widflx[ipnt-1]* 0.5 )
00406 {
00407 fprintf( ioQQQ, " ChckFill low test high fail\n" );
00408 lgFail = true;
00409 }
00410
00411
00412 energy = continuum.filbnd[i]*0.01 + continuum.filbnd[i+1]*0.99;
00413 ipnt = ipoint(energy);
00414
00415 if( energy < rfield.anu[ipnt-1] - rfield.widflx[ipnt-1]*0.5 )
00416 {
00417 fprintf( ioQQQ, " ChckFill high test low fail\n" );
00418 lgFail = true;
00419 }
00420
00421
00422
00423 else if( (energy > rfield.anu[ipnt-1] + rfield.widflx[ipnt-1]*0.5) &&
00424 ( energy > rfield.anu[ipnt] - rfield.widflx[ipnt]*0.5 ) )
00425 {
00426 fprintf( ioQQQ, " ChckFill high test high fail\n" );
00427 lgFail = true;
00428 }
00429 }
00430
00431 if( lgFail )
00432 {
00433 cdEXIT(EXIT_FAILURE);
00434 }
00435 return;
00436 }
00437
00438
00439 STATIC void rfield_opac_malloc(void)
00440 {
00441 long i;
00442
00443 DEBUG_ENTRY( "rfield_opac_malloc()" );
00444
00445
00446
00447 ++rfield.nupper;
00448
00449
00454
00455 rfield.fine_ener_lo = rfield.emm;
00456 rfield.fine_ener_hi = 1500.f;
00457
00458
00459
00460
00461
00462
00463
00464
00465
00466
00467
00468
00469
00470
00471
00474 double TeLowestFineOpacity = 1e4;
00475 rfield.fine_opac_velocity_width =
00476 (realnum)sqrt(2.*BOLTZMANN/ATOMIC_MASS_UNIT*TeLowestFineOpacity/
00477 dense.AtomicWeight[rfield.fine_opac_nelem] ) /
00478
00479
00480 rfield.fine_opac_nresolv;
00481
00482
00483 rfield.ipFineConVelShift = 0;
00484
00485
00486 rfield.fine_resol = rfield.fine_opac_velocity_width / SPEEDLIGHT;
00487
00488
00489 rfield.nfine_malloc = (long)(log10( rfield.fine_ener_hi / rfield.fine_ener_lo ) / log10( 1. + rfield.fine_resol ) );
00490 if( rfield.nfine_malloc <= 0 )
00491 TotalInsanity();
00492 rfield.nfine = rfield.nfine_malloc;
00493
00494
00495 rfield.fine_opac_zone = (realnum *)MALLOC(sizeof(realnum)*(unsigned)rfield.nfine_malloc );
00496 memset(rfield.fine_opac_zone , 0 , (unsigned long)rfield.nfine_malloc*sizeof(realnum) );
00497
00498
00499 rfield.fine_opt_depth = (realnum *)MALLOC(sizeof(realnum)*(unsigned)rfield.nfine_malloc );
00500 memset(rfield.fine_opt_depth , 0 , (unsigned long)rfield.nfine_malloc*sizeof(realnum) );
00501
00502 rfield.fine_anu = (realnum *)MALLOC(sizeof(realnum)*(unsigned)rfield.nfine_malloc );
00503
00504
00505 ASSERT( rfield.fine_ener_lo > 0. && rfield.fine_resol > 0 );
00506 for( i=0;i<rfield.nfine_malloc; ++i )
00507 {
00508 rfield.fine_anu[i] = rfield.fine_ener_lo * (realnum)pow( (1.+rfield.fine_resol), (i+1.) );
00509 }
00510
00511
00512
00513 rfield.line_count = (long *)MALLOC(sizeof(long)*(unsigned)NCELL );
00514 for( i=0; i<rfield.nupper; ++i)
00515 {
00516 rfield.line_count[i] = 0;
00517 }
00518 rfield.anu = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
00519 rfield.AnuOrg = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
00520 rfield.widflx = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
00521 rfield.anulog = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
00522 rfield.anusqr = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
00523 rfield.anu2 = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
00524 rfield.anu3 = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
00525 rfield.flux_beam_time = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
00526 rfield.flux_isotropic = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
00527 rfield.flux_beam_const = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
00528 rfield.flux_accum = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
00529 rfield.ExtinguishFactor = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
00530 rfield.convoc = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
00531 rfield.OccNumbBremsCont = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
00532 rfield.OccNumbIncidCont = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
00533 rfield.OccNumbDiffCont = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
00534 rfield.OccNumbContEmitOut = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
00535 rfield.ConInterOut = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
00536 rfield.SummedCon = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
00537 rfield.SummedDif = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
00538 rfield.SummedDifSave = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
00539 rfield.SummedOcc = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
00540 rfield.ConOTS_local_photons = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
00541 rfield.DiffuseEscape = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
00542 rfield.TotDiff2Pht = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
00543 rfield.ConOTS_local_OTS_rate = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
00544 rfield.otslin = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
00545 rfield.otscon = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
00546 rfield.outlin_noplot = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
00547 rfield.flux_beam_const_save = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
00548 rfield.flux_time_beam_save = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
00549 rfield.flux_isotropic_save = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
00550 rfield.trans_coef_zone = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
00551 rfield.trans_coef_total = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
00552 rfield.DiffuseLineEmission = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
00553 rfield.ipnt_coarse_2_fine = (long int*)MALLOC((size_t)(rfield.nupper*sizeof(long int)) );
00554
00555
00556 rfield.flux = (realnum**)MALLOC((size_t)(2*sizeof(realnum*)) );
00557 rfield.ConEmitReflec = (realnum**)MALLOC((size_t)(2*sizeof(realnum*)) );
00558 rfield.ConEmitOut = (realnum**)MALLOC((size_t)(2*sizeof(realnum*)) );
00559 rfield.ConRefIncid = (realnum**)MALLOC((size_t)(2*sizeof(realnum*)) );
00560 rfield.flux_total_incident = (realnum**)MALLOC((size_t)(2*sizeof(realnum*)) );
00561 rfield.reflin = (realnum**)MALLOC((size_t)(2*sizeof(realnum*)) );
00562 rfield.outlin = (realnum**)MALLOC((size_t)(2*sizeof(realnum*)) );
00563
00564 for( i=0; i<2; ++i )
00565 {
00566 rfield.flux[i] = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
00567 rfield.ConEmitReflec[i] = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
00568 rfield.ConEmitOut[i] = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
00569 rfield.ConRefIncid[i] = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
00570 rfield.flux_total_incident[i] = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
00571 rfield.reflin[i] = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
00572 rfield.outlin[i] = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
00573 }
00574
00575 memset(rfield.flux[1] , 0 , (unsigned long)rfield.nupper*sizeof(realnum) );
00576 memset(rfield.ConEmitReflec[1] , 0 , (unsigned long)rfield.nupper*sizeof(realnum) );
00577 memset(rfield.ConEmitOut[1] , 0 , (unsigned long)rfield.nupper*sizeof(realnum) );
00578 memset(rfield.ConRefIncid[1] , 0 , (unsigned long)rfield.nupper*sizeof(realnum) );
00579 memset(rfield.flux_total_incident[1] , 0 , (unsigned long)rfield.nupper*sizeof(realnum) );
00580 memset(rfield.reflin[1] , 0 , (unsigned long)rfield.nupper*sizeof(realnum) );
00581 memset(rfield.outlin[1] , 0 , (unsigned long)rfield.nupper*sizeof(realnum) );
00582
00583
00584
00585 rfield.gff = (realnum**)MALLOC((size_t)((LIMELM+1)*sizeof(realnum*)) );
00586 for( i = 1; i <= LIMELM; i++ )
00587 {
00588 rfield.gff[i] = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
00589 }
00590
00591 rfield.csigh = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
00592 rfield.csigc = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
00593
00594 rfield.comdn = (double*)MALLOC((size_t)(rfield.nupper*sizeof(double)) );
00595 rfield.comup = (double*)MALLOC((size_t)(rfield.nupper*sizeof(double)) );
00596 rfield.ContBoltz = (double*)MALLOC((size_t)(rfield.nupper*sizeof(double)) );
00597
00598
00599 rfield.otssav = (realnum**)MALLOC((size_t)(rfield.nupper*sizeof(realnum*)));
00600 for( i=0; i<rfield.nupper; ++i)
00601 {
00602 rfield.otssav[i] = (realnum*)MALLOC(2*sizeof(realnum));
00603 }
00604
00605
00606
00607 rfield.chLineLabel = (char**)MALLOC((size_t)(rfield.nupper*sizeof(char*)));
00608 rfield.chContLabel = (char**)MALLOC((size_t)(rfield.nupper*sizeof(char*)));
00609
00610
00611 for( i=0; i<rfield.nupper; ++i)
00612 {
00613 rfield.chLineLabel[i] = (char*)MALLOC(5*sizeof(char));
00614 rfield.chContLabel[i] = (char*)MALLOC(5*sizeof(char));
00615 }
00616
00617 opac.TauAbsFace = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
00618 memset( opac.TauAbsFace , 0 , rfield.nupper*sizeof(realnum) );
00619
00620 opac.TauScatFace = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
00621 opac.E2TauAbsFace = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
00622 opac.E2TauAbsTotal = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
00623 opac.TauAbsTotal = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
00624 opac.E2TauAbsOut = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
00625 opac.ExpmTau = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
00626 opac.tmn = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
00627
00628 opac.opacity_abs = (double*)MALLOC((size_t)(rfield.nupper*sizeof(double)) );
00629 opac.opacity_abs_savzon1 = (double*)MALLOC((size_t)(rfield.nupper*sizeof(double)) );
00630 opac.OldOpacSave = (double*)MALLOC((size_t)(rfield.nupper*sizeof(double)) );
00631 opac.opacity_sct = (double*)MALLOC((size_t)(rfield.nupper*sizeof(double)) );
00632 opac.albedo = (double*)MALLOC((size_t)(rfield.nupper*sizeof(double)) );
00633 opac.opacity_sct_savzon1 = (double*)MALLOC((size_t)(rfield.nupper*sizeof(double)) );
00634 opac.OpacStatic = (double*)MALLOC((size_t)(rfield.nupper*sizeof(double)) );
00635 opac.FreeFreeOpacity = (double*)MALLOC((size_t)(rfield.nupper*sizeof(double)) );
00636 opac.ExpZone = (double*)MALLOC((size_t)(rfield.nupper*sizeof(double)) );
00637
00638 opac.TauAbsGeo = (realnum**)MALLOC((size_t)(2*sizeof(realnum *)) );
00639 opac.TauScatGeo = (realnum**)MALLOC((size_t)(2*sizeof(realnum *)) );
00640 opac.TauTotalGeo = (realnum**)MALLOC((size_t)(2*sizeof(realnum *)) );
00641
00642 for( i=0; i<2; ++i)
00643 {
00644 opac.TauAbsGeo[i] = (realnum*)MALLOC(rfield.nupper*sizeof(realnum));
00645 opac.TauScatGeo[i] = (realnum*)MALLOC(rfield.nupper*sizeof(realnum));
00646 opac.TauTotalGeo[i] = (realnum*)MALLOC(rfield.nupper*sizeof(realnum));
00647 }
00648
00649
00650 --rfield.nupper;
00651
00652
00653 lgRfieldMalloced = true;
00654 return;
00655 }
00656
00657
00658
00659 STATIC void read_continuum_mesh( void )
00660 {
00661 FILE *ioDATA;
00662 char chLine[INPUT_LINE_LENGTH];
00663 long i;
00664 bool lgEOL;
00665 long i1 , i2 , i3;
00666
00667 DEBUG_ENTRY( "read_continuum_mesh()" );
00668
00669 if( trace.lgTrace )
00670 fprintf( ioQQQ," read_continuum_mesh opening continuum_mesh.ini:");
00671
00672 ioDATA = open_data( "continuum_mesh.ini", "r" );
00673
00674
00675 if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
00676 {
00677 fprintf( ioQQQ, " read_continuum_mesh could not read first line of continuum_mesh.ini.\n");
00678 cdEXIT(EXIT_FAILURE);
00679 }
00680
00681
00682 continuum.nStoredBands = 0;
00683 while( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) != NULL )
00684 {
00685
00686
00687 if( chLine[0] != '#')
00688 ++continuum.nStoredBands;
00689 }
00690
00691
00692
00693 continuum.filbnd =
00694 ((realnum *)MALLOC( (size_t)(continuum.nStoredBands+1)*sizeof(realnum )));
00695 continuum.fildel =
00696 ((realnum *)MALLOC( (size_t)(continuum.nStoredBands+1)*sizeof(realnum )));
00697 continuum.filres =
00698 ((realnum *)MALLOC( (size_t)(continuum.nStoredBands+1)*sizeof(realnum )));
00699 continuum.ifill0 =
00700 ((long *)MALLOC( (size_t)(continuum.nStoredBands+1)*sizeof(long )));
00701 continuum.StoredEnergy =
00702 ((double *)MALLOC( (size_t)(continuum.nStoredBands+1)*sizeof(double )));
00703 continuum.StoredResolution =
00704 ((double *)MALLOC( (size_t)(continuum.nStoredBands+1)*sizeof(double )));
00705
00706
00707 if( fseek( ioDATA , 0 , SEEK_SET ) != 0 )
00708 {
00709 fprintf( ioQQQ, " read_continuum_mesh could not rewind continuum_mesh.ini.\n");
00710 cdEXIT(EXIT_FAILURE);
00711 }
00712
00713
00714 if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
00715 {
00716 fprintf( ioQQQ, " read_continuum_mesh could not read first line of continuum_mesh.ini.\n");
00717 cdEXIT(EXIT_FAILURE);
00718 }
00719
00720 i = 1;
00721
00722 i1 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
00723 i2 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
00724 i3 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
00725
00726 bool lgResPower;
00727
00728
00729 if( i1 == 1 && i2 == 9 && i3 == 29 )
00730
00731
00732 lgResPower = false;
00733 else if( i1 == 10 && i2 == 8 && i3 == 8 )
00734
00735
00736 lgResPower = true;
00737 else
00738 {
00739 fprintf( ioQQQ,
00740 " read_continuum_mesh: the version of continuum_mesh.ini is not supported.\n" );
00741 fprintf( ioQQQ,
00742 " I found version number %li %li %li.\n" ,
00743 i1 , i2 , i3 );
00744 fprintf( ioQQQ, "Here is the line image:\n==%s==\n", chLine );
00745 cdEXIT(EXIT_FAILURE);
00746 }
00747
00748
00749
00750 continuum.nStoredBands = 0;
00751 while( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) != NULL )
00752 {
00753
00754 if( chLine[0] != '#')
00755 {
00756 i = 1;
00757 continuum.StoredEnergy[continuum.nStoredBands] = FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
00758 continuum.StoredResolution[continuum.nStoredBands] = FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
00759
00760
00761
00762 if( continuum.StoredEnergy[continuum.nStoredBands]<0. ||
00763 continuum.StoredResolution[continuum.nStoredBands]<=0. )
00764 {
00765 fprintf(ioQQQ, "DISASTER PROBLEM continuum_mesh.ini has a non-positive number.\n");
00766 cdEXIT(EXIT_FAILURE);
00767 }
00768
00769
00770 if( lgResPower )
00771 continuum.StoredResolution[continuum.nStoredBands] = 1./
00772 continuum.StoredResolution[continuum.nStoredBands];
00773
00774
00775 continuum.StoredResolution[continuum.nStoredBands] *= continuum.ResolutionScaleFactor;
00776
00777 ++continuum.nStoredBands;
00778 }
00779 }
00780
00781 fclose( ioDATA );
00782
00783
00784 for( i=1; i<continuum.nStoredBands-1; ++i )
00785 {
00786 if( continuum.StoredEnergy[i-1]>=continuum.StoredEnergy[i] )
00787 {
00788 fprintf( ioQQQ,
00789 " read_continuum_mesh: The continuum definition array energies must be in increasing order.\n" );
00790 cdEXIT(EXIT_FAILURE);
00791 }
00792 }
00793 if( continuum.StoredEnergy[continuum.nStoredBands-1]!=0 )
00794 {
00795 fprintf( ioQQQ,
00796 " read_continuum_mesh: The last continuum array energies must be zero.\n" );
00797 cdEXIT(EXIT_FAILURE);
00798 }
00799 return;
00800 }
00801
00802
00803 void rfield_opac_zero(
00804
00805 long lo ,
00806
00807 long ihi )
00808 {
00809 long int i;
00810
00811
00812
00813
00814 if( lgRfieldMalloced )
00815 {
00816 unsigned long n=(unsigned long)(ihi-lo+1);
00817 memset(&rfield.OccNumbDiffCont[lo] , 0 , n*sizeof(realnum) );
00818 memset(&rfield.OccNumbContEmitOut[lo] , 0 , n*sizeof(realnum) );
00819 memset(&rfield.ContBoltz[lo] , 0 , n*sizeof(double) );
00820
00821
00822
00823 memset(&rfield.ConEmitReflec[0][lo] , 0 , n*sizeof(realnum) );
00824 memset(&rfield.ConEmitOut[0][lo] , 0 , n*sizeof(realnum) );
00825 memset(&rfield.reflin[0][lo] , 0 , n*sizeof(realnum) );
00826 memset(&rfield.ConRefIncid[0][lo] , 0 , n*sizeof(realnum) );
00827 memset(&rfield.SummedCon[lo] , 0 , n*sizeof(realnum) );
00828 memset(&rfield.OccNumbBremsCont[lo] , 0 , n*sizeof(realnum) );
00829 memset(&rfield.convoc[lo] , 0 , n*sizeof(realnum) );
00830 memset(&rfield.flux[0][lo] , 0 , n*sizeof(realnum) );
00831 memset(&rfield.flux_total_incident[0][lo] , 0 , n*sizeof(realnum) );
00832 memset(&rfield.flux_beam_const_save[lo] , 0 , n*sizeof(realnum) );
00833 memset(&rfield.flux_time_beam_save[lo] , 0 , n*sizeof(realnum) );
00834 memset(&rfield.flux_isotropic_save[lo] , 0 , n*sizeof(realnum) );
00835 memset(&rfield.SummedOcc[lo] , 0 , n*sizeof(realnum) );
00836 memset(&rfield.SummedDif[lo] , 0 , n*sizeof(realnum) );
00837 memset(&rfield.flux_accum[lo] , 0 , n*sizeof(realnum) );
00838 memset(&rfield.otslin[lo] , 0 , n*sizeof(realnum) );
00839 memset(&rfield.otscon[lo] , 0 , n*sizeof(realnum) );
00840 memset(&rfield.ConInterOut[lo] , 0 , n*sizeof(realnum) );
00841 memset(&rfield.outlin[0][lo] , 0 , n*sizeof(realnum) );
00842 memset(&rfield.outlin_noplot[lo] , 0 , n*sizeof(realnum) );
00843 memset(&rfield.ConOTS_local_OTS_rate[lo], 0 , n*sizeof(realnum) );
00844 memset(&rfield.ConOTS_local_photons[lo] , 0 , n*sizeof(realnum) );
00845 memset(&opac.OldOpacSave[lo] , 0 , n*sizeof(double) );
00846 memset(&opac.opacity_abs[lo] , 0 , n*sizeof(double) );
00847 memset(&opac.opacity_sct[lo] , 0 , n*sizeof(double) );
00848 memset(&opac.albedo[lo] , 0 , n*sizeof(double) );
00849 memset(&opac.FreeFreeOpacity[lo] , 0 , n*sizeof(double) );
00850
00851
00852 memset( &opac.E2TauAbsTotal[lo] , 0 , n*sizeof(realnum) );
00853 memset( &opac.E2TauAbsOut[lo] , 0 , n*sizeof(realnum) );
00854 memset( &opac.TauAbsTotal[lo] , 0 , n*sizeof(realnum) );
00855
00856 for( i=lo; i <= ihi; i++ )
00857 {
00858 opac.TauTotalGeo[0][i] = opac.taumin;
00859 opac.TauAbsGeo[0][i] = opac.taumin;
00860 opac.TauScatGeo[0][i] = opac.taumin;
00861 opac.tmn[i] = 1.;
00862 opac.ExpZone[i] = 1.;
00863 opac.E2TauAbsFace[i] = 1.;
00864 opac.ExpmTau[i] = 1.;
00865 opac.OpacStatic[i] = 1.;
00866 }
00867
00868 memset(rfield.fine_opac_zone , 0 , (unsigned long)rfield.nfine_malloc*sizeof(realnum) );
00869
00870 memset(rfield.fine_opt_depth , 0 , (unsigned long)rfield.nfine_malloc*sizeof(realnum) );
00871 }
00872 return;
00873 }