cloudy
trunk
|
00001 /* This file is part of Cloudy and is copyright (C)1978-2008 by Gary J. Ferland and 00002 * others. For conditions of distribution and use see copyright notice in license.txt */ 00003 /*ContCreateMesh calls fill to set up continuum energy mesh if first call, 00004 * otherwise reset to original mesh */ 00005 /*fill define the continuum energy grid over a specified range */ 00006 /*ChckFill perform sanity check confirming that the energy array has been properly filled */ 00007 /*rfield_opac_malloc MALLOC space for opacity arrays */ 00008 /*read_continuum_mesh read the continuum definition from the file continuum_mesh.ini */ 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 /* read the continuum definition from the file continuum_mesh.ini */ 00022 STATIC void read_continuum_mesh( void ); 00023 00024 /*fill define the continuum energy grid over a specified range */ 00025 STATIC void fill(double fenlo, 00026 double fenhi, 00027 double resolv, 00028 long int *n0, 00029 long int *ipnt, 00030 /* this says only count, do not fill */ 00031 bool lgCount ); 00032 00033 /*rfield_opac_malloc MALLOC space for opacity arrays */ 00034 STATIC void rfield_opac_malloc(void); 00035 00036 /*ChckFill perform sanity check confirming that the energy array has been properly filled */ 00037 STATIC void ChckFill(void); 00038 00039 void ContCreateMesh(void) 00040 { 00041 long int 00042 i, 00043 ipnt, 00044 n0; 00045 00046 /* flag to say whether pointers have ever been evaluated */ 00047 static bool lgPntEval = false; 00048 00049 DEBUG_ENTRY( "ContCreateMesh()" ); 00050 00051 /* lgPntEval is local static variable defined false when defined. 00052 * it is set true below, so that pointers only created one time in the 00053 * history of this coreload. */ 00054 if( lgPntEval ) 00055 { 00056 if( trace.lgTrace ) 00057 { 00058 fprintf( ioQQQ, " ContCreateMesh called, not evaluating.\n" ); 00059 } 00060 /* now save current form of energy array */ 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 /* set size of arrays that continuum iteration information */ 00078 00079 /* read in the continuum mesh resolution definition */ 00080 /* >>chng 01 sep 29, add external file "continuum_mesh.ini" with fill parameters */ 00081 read_continuum_mesh(); 00082 00083 /* fill in continuum with freq points 00084 * arg are range pointer, 2 energy limits, resolution 00085 * first argument is lower energy of the continuum range to be defined 00086 * second argument is upper range; both are in Rydbergs 00087 * third number is the relative energy resolution, dnu/nu, 00088 * for that range of the continuum 00089 * last two numbers are internal book keeping 00090 * N0 is number of energy cells used so far 00091 * IPNT is a counter for the number of fills used 00092 * 00093 * if this is changed, then also change warning in GetTable about using 00094 * transmitted continuum - it says version number where continuum changed 00095 * */ 00096 n0 = 1; 00097 ipnt = 0; 00098 /* this is number of ranges that will be introduced below*/ 00099 continuum.nrange = 0; 00100 00101 /* ================================================================ */ 00102 /* NB - this block must agree exactly with the one that follows */ 00103 n0 = 1; 00104 ipnt = 0; 00105 /* this is number of ranges that will be introduced below*/ 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 /* at this point debugger shows that anu and widflx are defined 00120 * up through n0-2 - due to c offset -1 from Fortran! */ 00121 rfield.nupper = n0 - 1; 00122 /* there must be a cell above nflux for us to pass unity through the vol integrator */ 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 /*>>chng 04 oct 10, from nflux = nupper to nupper-1 since vectors must malloc to nupper, but 00132 * will address [nflux] for unit continuum test */ 00133 rfield.nflux = rfield.nupper-1; 00134 00135 /* allocate space for continuum arrays within rfield.h and opacity arrays in opacity.h 00136 * sets lgRfieldMalloced true */ 00137 rfield_opac_malloc(); 00138 00139 /* geometry.nend_max is largest number of zones needed on any iteration, 00140 * will use it to malloc arrays that save source function as function of zone */ 00141 /* now change all limits, for all iterations, to this value */ 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 /* nend_max+1 because search phase is zone 0, first zone at illumin face is 1 */ 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 /* this is number of ranges that will be introduced below*/ 00159 continuum.nrange = 0; 00160 00161 /* the default array values are set in continuum_mesh.ini */ 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 /* fill in the false highest cell used for unit verification */ 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 /* there must be a cell above nflux for us to pass unity through the vol integrator 00180 * as a sanity check. assert that this is true so we will crash if ever changed 00181 ASSERT( rfield.nupper +1 <= rfield.nupper );*/ 00182 00183 /* this is done here when the space is first allocated, 00184 * then done on every subsequent initialization in zero.c */ 00185 rfield_opac_zero( 0 , rfield.nupper ); 00186 00187 /* this is a sanity check for results produced above by fill */ 00188 ChckFill(); 00189 00190 /* now fix widflx array so that it is correct */ 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 /* now save current form of array, and define some quantities related to it */ 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 /* following are Compton exchange factors from Tarter */ 00206 /* this code also appears in highen, but coef needed before that routine called. */ 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 /* >>chng 05 feb 28, add transmission and mapping coef */ 00214 /* map these coarse continua into fine continuum grid */ 00215 if( rfield.anu[i] < rfield.fine_ener_lo || rfield.anu[i]>rfield.fine_ener_hi ) 00216 { 00217 /* 0 (false) says not defined */ 00218 rfield.ipnt_coarse_2_fine[i] = 0; 00219 } 00220 else 00221 { 00222 if( ipnt==0 ) 00223 { 00224 /* this is the first one that maps onto the fine array */ 00225 rfield.ipnt_coarse_2_fine[i] = 0; 00226 ipnt = 1; 00227 } 00228 else 00229 { 00230 /* find first fine frequency that is greater than this coarse value */ 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 /*fprintf(ioQQQ," coarse %li nu= %.3e points to fine %li nu=%.3e\n", 00239 i, rfield.anu[i] , rfield.ipnt_coarse_2_fine[i] , rfield.fine_anu[rfield.ipnt_coarse_2_fine[i]] );*/ 00240 rfield.trans_coef_zone[i] = 1.f; 00241 rfield.trans_coef_total[i] = 1.f; 00242 } 00243 return; 00244 } 00245 00246 /*fill define the continuum energy grid over a specified range, called by ContCreateMesh */ 00247 STATIC void fill( 00248 /* lower bounds to this energy range */ 00249 double fenlo, 00250 /* upper bounds to this continuum range */ 00251 double fenhi, 00252 /* relative energy resolution */ 00253 double resolv, 00254 /* starting index within frequency grid */ 00255 long int *n0, 00256 /* which energy band this is */ 00257 long int *ipnt, 00258 /* this says only count, do not fill */ 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 /* this is the number of cells needed to fill the array with numbers at the requested resolution */ 00271 nbin = (long int)(log(10.)*log10(fenhi/fenlo)/resolv + 1); 00272 00273 if( lgCount ) 00274 { 00275 /* true means only count number of cells, don't do anything */ 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 /* this is the number of cells needed to fill the array with numbers 00295 nbin = (long int)(log(10.)*log10(fenhi/fenlo)/resolv + 1);*/ 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 /* nrange is number of ranges */ 00348 *ipnt += 1; 00349 continuum.nrange = MAX2(continuum.nrange,*ipnt); 00350 return; 00351 } 00352 00353 /*ChckFill perform sanity check confirming that the energy array has been properly filled */ 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 /* test middle of energy bound */ 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 /* >>chng 02 jul 16, add second test - when "set resol 10" used, 00379 * very large values of cell width, combined with fact that cells 00380 * are log increasing, causes problem. */ 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 /* test near low bound */ 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 /* test near high bound */ 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 /* >>chng 02 jul 16, add second test - when "set resol 10" used, 00413 * very large values of cell width, combined with fact that cells 00414 * are log increasing, causes problem. */ 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 /* MALLOC arrays within rfield */ 00431 STATIC void rfield_opac_malloc(void) 00432 { 00433 long i; 00434 00435 DEBUG_ENTRY( "rfield_opac_malloc()" ); 00436 00437 /* allocate one more than we use for the unit integration, 00438 * will back up at end of routine */ 00439 ++rfield.nupper; 00440 00441 /* >>chng 03 feb 12, add fine mesh fine grid fine opacity array to keep track of line overlap */ 00446 /* frequency range in Rydberg needed for all resonance lines */ 00447 rfield.fine_ener_lo = 0.05f; 00448 rfield.fine_ener_hi = 1500.f; 00449 00450 /* set resolution of fine continuum mesh. 00451 * rfield.fine_opac_velocity_width is width per cell, cm/s 00452 * choose width so that most massive species (usually Fe) is well resolved 00453 * 00454 * rfield.fine_opac_nelem is the most massive (hence sharpest line) 00455 * we will worry about. By default this is irion but can be changed 00456 * with SET FINE CONTINUUM command 00457 * 00458 * TeLowestFineOpacity of 1e4 K is temperature were line width is 00459 * evaluated. Tests were done using the stop temperature in its place 00460 * Te below 1e4 K made fine opacity grid huge 00461 * do not let temp get higher than 1e4 either - code run with stop temp 10 set 00462 * stop temp of 1e10K and assert thrown at line 204 of cont_createpointers.c 00463 * simply use 1e4 K as a characteristic temperature */ 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 /* we want fine_opac_nresolv continuum elements across this line 00471 * default is 1, changed with SET FINE CONTINUUM command */ 00472 rfield.fine_opac_nresolv; 00473 00474 /* we are at first zone so velocity shift is zero */ 00475 rfield.ipFineConVelShift = 0; 00476 00477 /* dimensionless resolution, dE/E, this is used in ipoint to get offset in find mesh */ 00478 rfield.fine_resol = rfield.fine_opac_velocity_width / SPEEDLIGHT; 00479 00480 /* the number of cells needed */ 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 /* this is the fine opacity array to ghost the main low-resolution array */ 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 /* this is the fine total optical array to ghost the main low-resolution array */ 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 /* now fill in energy array */ 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 /* done with fine array */ 00503 00504 /* used to count number of lines per cell */ 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 /* chng 02 may 16, by Ryan...added array for gaunt factors for ALL charges, malloc here. */ 00554 /* First index is EFFECTIVE CHARGE MINUS ONE! */ 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 /*realnum rfield.otssav[NC_ELL][2];*/ 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 /* char rfield.chLineLabel[NLINES][5];*/ 00578 rfield.chLineLabel = (char**)MALLOC((size_t)(rfield.nupper*sizeof(char*))); 00579 rfield.chContLabel = (char**)MALLOC((size_t)(rfield.nupper*sizeof(char*))); 00580 00581 /* now allocate all the labels for each of the above lines */ 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 /* fix allocate trick for one more than we use for the unit integration */ 00621 --rfield.nupper; 00622 00623 /* say that space exists */ 00624 lgRfieldMalloced = true; 00625 return; 00626 } 00627 00628 00629 /* read the continuum definition from the file continuum_mesh.ini */ 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 /* first line is a version number and does not count */ 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 /* count how many lines are in the file, ignoring all lines 00652 * starting with '#' */ 00653 continuum.nStoredBands = 0; 00654 while( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) != NULL ) 00655 { 00656 /* we want to count the lines that do not start with # 00657 * since these contain data */ 00658 if( chLine[0] != '#') 00659 ++continuum.nStoredBands; 00660 } 00661 00662 /* we now have number of lines containing pairs of bounds, 00663 * allocate space for the arrays we will need */ 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 /* now rewind the file so we can read it a second time*/ 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 /* check that magic number is ok */ 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 /* level 1 magic number */ 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 /* the following is the set of numbers that appear at the start of continuum_mesh.ini 01 08 10 */ 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 /* this starts at 1 not 0 since zero is reserved for the 00710 * dummy line */ 00711 continuum.nStoredBands = 0; 00712 while( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) != NULL ) 00713 { 00714 /* only look at lines without '#' in first col */ 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 /* option to enter numbers as logs if less than zero */ 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 /* this is option to rescale resolution with set resolution command */ 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 /* now verify continuum grid is ok - first are all values but the last positive? */ 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 }