cloudy trunk
Loading...
Searching...
No Matches
cont_createmesh.cpp
Go to the documentation of this file.
1/* This file is part of Cloudy and is copyright (C)1978-2013 by Gary J. Ferland and
2 * others. For conditions of distribution and use see copyright notice in license.txt */
3/*ContCreateMesh calls fill to set up continuum energy mesh if first call,
4 * otherwise reset to original mesh */
5/*fill define the continuum energy grid over a specified range */
6/*ChckFill perform sanity check confirming that the energy array has been properly filled */
7/*rfield_opac_malloc MALLOC space for opacity arrays */
8/*read_continuum_mesh read the continuum definition from the file continuum_mesh.ini */
9#include "cddefines.h"
10#include "rfield.h"
11#include "iterations.h"
12#include "physconst.h"
13#include "dense.h"
14#include "trace.h"
15#include "opacity.h"
16#include "ipoint.h"
17#include "geometry.h"
18#include "continuum.h"
19
20/* read the continuum definition from the file continuum_mesh.ini */
21STATIC void read_continuum_mesh( void );
22
23/*fill define the continuum energy grid over a specified range */
24STATIC void fill(double fenlo,
25 double fenhi,
26 double resolv,
27 long int *n0,
28 long int *ipnt,
29 /* this says only count, do not fill */
30 bool lgCount );
31
32/*rfield_opac_malloc MALLOC space for opacity arrays */
33STATIC void rfield_opac_malloc(void);
34
35/*ChckFill perform sanity check confirming that the energy array has been properly filled */
36STATIC void ChckFill(void);
37
39{
40 long int
41 i,
42 ipnt,
43 n0;
44
45 /* flag to say whether pointers have ever been evaluated */
46 static bool lgPntEval = false;
47
48 DEBUG_ENTRY( "ContCreateMesh()" );
49
50 /* lgPntEval is local static variable defined false when defined.
51 * it is set true below, so that pointers only created one time in the
52 * history of this coreload. */
53 if( lgPntEval )
54 {
55 if( trace.lgTrace )
56 {
57 fprintf( ioQQQ, " ContCreateMesh called, not evaluating.\n" );
58 }
59 /* now save current form of energy array */
60 for( i=0; i < rfield.nupper; i++ )
61 {
62 rfield.anu[i] = rfield.AnuOrg[i];
63 rfield.anu2[i] = rfield.anu[i]*rfield.anu[i];
64 }
65 memset( opac.TauAbsFace , 0 , rfield.nupper*sizeof(realnum) );
66 return;
67 }
68 else
69 {
70 if( trace.lgTrace )
71 {
72 fprintf( ioQQQ, " ContCreateMesh called first time.\n" );
73 }
74 lgPntEval = true;
75 }
76
77 /* set size of arrays that continuum iteration information */
78
79 /* read in the continuum mesh resolution definition */
80 /* >>chng 01 sep 29, add external file "continuum_mesh.ini" with fill parameters */
82
83 /* fill in continuum with freq points
84 * arg are range pointer, 2 energy limits, resolution
85 * first argument is lower energy of the continuum range to be defined
86 * second argument is upper range; both are in Rydbergs
87 * third number is the relative energy resolution, dnu/nu,
88 * for that range of the continuum
89 * last two numbers are internal book keeping
90 * N0 is number of energy cells used so far
91 * IPNT is a counter for the number of fills used
92 *
93 * if this is changed, then also change warning in GetTable about using
94 * transmitted continuum - it says version number where continuum changed
95 * */
96 n0 = 1;
97 ipnt = 0;
98 /* this is number of ranges that will be introduced below*/
99 continuum.nrange = 0;
100
101 /* ================================================================ */
102 /* NB - this block must agree exactly with the one that follows */
103 n0 = 1;
104 ipnt = 0;
105 /* this is number of ranges that will be introduced below*/
106 continuum.nrange = 0;
107 continuum.StoredEnergy[continuum.nStoredBands-1] = rfield.egamry;
108
109 fill(rfield.emm, continuum.StoredEnergy[0] , continuum.StoredResolution[0],&n0,&ipnt,true );
110 for(i=1; i<continuum.nStoredBands; ++i )
111 {
112 fill(continuum.StoredEnergy[i-1] ,
113 continuum.StoredEnergy[i],
114 continuum.StoredResolution[i],
115 &n0,&ipnt,true);
116 }
117 /* ================================================================ */
118
119 /* at this point debugger shows that anu and widflx are defined
120 * up through n0-2 - due to c offset -1 from Fortran! */
121 rfield.nupper = n0 - 1;
122 /* there must be a cell above nflux for us to pass unity through the vol integrator */
123 if( rfield.nupper >= NCELL )
124 {
125 fprintf(ioQQQ," Currently the arrays that hold interpolated tables can only hold %i points.\n",NCELL);
126 fprintf(ioQQQ," This continuum mesh really needs to have %li points.\n",rfield.nupper);
127 fprintf(ioQQQ," Please increase the value of NCELL in rfield.h and recompile.\n Sorry.");
129 }
130
131 /*>>chng 04 oct 10, from nflux = nupper to nupper-1 since vectors must malloc to nupper, but
132 * will address [nflux] for unit continuum test */
133 rfield.nflux = rfield.nupper-1;
134
135 /* allocate space for continuum arrays within rfield.h and opacity arrays in opacity.h
136 * sets lgRfieldMalloced true */
138
139 /* geometry.nend_max is largest number of zones needed on any iteration,
140 * will use it to malloc arrays that save source function as function of zone */
141 /* now change all limits, for all iterations, to this value */
142 geometry.nend_max = geometry.nend[0];
143 for( i=1; i < iterations.iter_malloc; i++ )
144 {
145 geometry.nend_max = MAX2( geometry.nend_max , geometry.nend[i] );
146 }
147 /* nend_max+1 because search phase is zone 0, first zone at illumin face is 1 */
148 rfield.ConEmitLocal = (realnum**)MALLOC( (size_t)(geometry.nend_max+1)*sizeof(realnum *) );
149 rfield.ConSourceFcnLocal = (realnum**)MALLOC( (size_t)(geometry.nend_max+1)*sizeof(realnum *) );
150 for( i=0;i<(geometry.nend_max+1); ++i )
151 {
152 rfield.ConEmitLocal[i] = (realnum*)MALLOC( (size_t)rfield.nupper*sizeof(realnum) );
153 rfield.ConSourceFcnLocal[i] = (realnum*)MALLOC( (size_t)rfield.nupper*sizeof(realnum) );
154 }
155 for( i=0;i<(geometry.nend_max+1); ++i )
156 {
157 for( long j=0; j<rfield.nupper; ++j)
158 {
159 rfield.ConSourceFcnLocal[i][j] = 1.;
160 }
161 }
162
163 /* ================================================================ */
164 n0 = 1;
165 ipnt = 0;
166
167 /* this is number of ranges that will be introduced below*/
168 continuum.nrange = 0;
169
170 /* the default array values are set in continuum_mesh.ini */
171 fill(rfield.emm, continuum.StoredEnergy[0] , continuum.StoredResolution[0] ,
172 &n0,&ipnt,false);
173 for(i=1; i<continuum.nStoredBands; ++i )
174 {
175 fill(continuum.StoredEnergy[i-1] ,
176 continuum.StoredEnergy[i],
177 continuum.StoredResolution[i],
178 &n0,&ipnt,false);
179 }
180
181 /* ================================================================ */
182
183 /* fill in the false highest cell used for unit verification */
184 rfield.widflx[rfield.nupper] = rfield.widflx[rfield.nupper-1];
185 rfield.anu[rfield.nupper] = rfield.anu[rfield.nupper-1] +
186 rfield.widflx[rfield.nupper];
187
188 /* there must be a cell above nflux for us to pass unity through the vol integrator
189 * as a sanity check. assert that this is true so we will crash if ever changed
190 ASSERT( rfield.nupper +1 <= rfield.nupper );*/
191
192 /* this is done here when the space is first allocated,
193 * then done on every subsequent initialization in zero.c */
194 rfield_opac_zero( 0 , rfield.nupper );
195
196 /* this is a sanity check for results produced above by fill */
197 ChckFill();
198
199 /* now fix widflx array so that it is correct */
200 for( i=1; i<rfield.nupper-1; ++i )
201 {
202 rfield.widflx[i] = ((rfield.anu[i+1] - rfield.anu[i]) + (rfield.anu[i] -
203 rfield.anu[i-1]))/2.f;
204 }
205
206 ipnt = 0;
207 /* now save current form of array, and define some quantities related to it */
208 for( i=0; i < rfield.nupper; i++ )
209 {
210 double alf , bet;
211
212 rfield.AnuOrg[i] = rfield.anu[i];
213 rfield.anusqr[i] = (realnum)sqrt(rfield.AnuOrg[i]);
214 /* following are Compton exchange factors from Tarter */
215 /* this code also appears in highen, but coef needed before that routine called. */
216 alf = 1./(1. + rfield.anu[i]*(1.1792e-4 + 7.084e-10*rfield.anu[i]));
217 bet = 1. - alf*rfield.anu[i]*(1.1792e-4 + 2.*7.084e-10*rfield.anu[i])/4.;
218 rfield.csigh[i] = (realnum)(alf*rfield.anu[i]*rfield.anu[i]*3.858e-25);
219 rfield.csigc[i] = (realnum)(alf*bet*rfield.anu[i]*3.858e-25);
220 rfield.anu2[i] = rfield.anu[i]*rfield.anu[i];
221
222 /* >>chng 05 feb 28, add transmission and mapping coef */
223 /* map these coarse continua into fine continuum grid */
224 if( rfield.anu[i] < rfield.fine_ener_lo || rfield.anu[i]>rfield.fine_ener_hi )
225 {
226 /* 0 (false) says not defined */
227 rfield.ipnt_coarse_2_fine[i] = 0;
228 }
229 else
230 {
231 if( ipnt==0 )
232 {
233 /* this is the first one that maps onto the fine array */
234 rfield.ipnt_coarse_2_fine[i] = 0;
235 ipnt = 1;
236 }
237 else
238 {
239 /* find first fine frequency that is greater than this coarse value */
240 while (ipnt < rfield.nfine_malloc && rfield.fine_anu[ipnt]<rfield.anu[i] )
241 {
242 ++ipnt;
243 }
244 rfield.ipnt_coarse_2_fine[i] = ipnt;
245 }
246 }
247 /*fprintf(ioQQQ," coarse %li nu= %.3e points to fine %li nu=%.3e\n",
248 i, rfield.anu[i] , rfield.ipnt_coarse_2_fine[i] , rfield.fine_anu[rfield.ipnt_coarse_2_fine[i]] );*/
249 }
250 rfield.resetCoarseTransCoef();
251 return;
252}
253
254/*fill define the continuum energy grid over a specified range, called by ContCreateMesh */
256 /* lower bounds to this energy range */
257 double fenlo,
258 /* upper bounds to this continuum range */
259 double fenhi,
260 /* relative energy resolution */
261 double resolv,
262 /* starting index within frequency grid */
263 long int *n0,
264 /* which energy band this is */
265 long int *ipnt,
266 /* this says only count, do not fill */
267 bool lgCount )
268{
269 long int i,
270 nbin;
271 realnum widtot;
272 double aaa , bbb;
273
274 DEBUG_ENTRY( "fill()" );
275
276 ASSERT( fenlo>0. && fenhi>0. && resolv>0. );
277
278 /* this is the number of cells needed to fill the array with numbers at the requested resolution */
279 nbin = (long int)(log(10.)*log10(fenhi/fenlo)/resolv + 1);
280
281 if( lgCount )
282 {
283 /* true means only count number of cells, don't do anything */
284 *n0 += nbin;
285 return;
286 }
287
288 if( *ipnt > 0 && fabs(1.-fenlo/continuum.filbnd[*ipnt]) > 1e-4 )
289 {
290 fprintf( ioQQQ, " FILL improper bounds.\n" );
291 fprintf( ioQQQ, " ipnt=%3ld fenlo=%11.4e filbnd(ipnt)=%11.4e\n",
292 *ipnt, fenlo, continuum.filbnd[*ipnt] );
294 }
295
296 ASSERT( *ipnt < continuum.nStoredBands );
297
298 continuum.ifill0[*ipnt] = *n0 - 1;
299 continuum.filbnd[*ipnt] = (realnum)fenlo;
300 continuum.filbnd[*ipnt+1] = (realnum)fenhi;
301
302 /* this is the number of cells needed to fill the array with numbers
303 nbin = (long int)(log(10.)*log10(fenhi/fenlo)/resolv + 1);*/
304 continuum.fildel[*ipnt] = (realnum)(log10(fenhi/fenlo)/nbin);
305
306 if( continuum.fildel[*ipnt] < 0.01 )
307 {
308 continuum.filres[*ipnt] = (realnum)(log(10.)*continuum.fildel[*ipnt]);
309 }
310 else
311 {
312 continuum.filres[*ipnt] = (realnum)((pow(10.,2.*continuum.fildel[*ipnt]) - 1.)/2./
313 pow((realnum)10.f,continuum.fildel[*ipnt]));
314 }
315
316 if( (*n0 + nbin-2) > rfield.nupper )
317 {
318 fprintf( ioQQQ, " Fill would need %ld cells to get to an energy of %.3e\n",
319 *n0 + nbin, fenhi );
320 fprintf( ioQQQ, " This is a major logical error in fill.\n");
321 ShowMe();
323 }
324
325 widtot = 0.;
326 for( i=0; i < nbin; i++ )
327 {
328 bbb = continuum.fildel[*ipnt]*((realnum)(i) + 0.5);
329 aaa = pow( 10. , bbb );
330
331 rfield.anu[i+continuum.ifill0[*ipnt]] = (realnum)(fenlo*aaa);
332
333 rfield.widflx[i+continuum.ifill0[*ipnt]] = rfield.anu[i+continuum.ifill0[*ipnt]]*
334 continuum.filres[*ipnt];
335
336 widtot += rfield.widflx[i+continuum.ifill0[*ipnt]];
337 }
338
339 *n0 += nbin;
340 if( trace.lgTrace && (trace.lgConBug || trace.lgPtrace) )
341 {
342 fprintf( ioQQQ,
343 " FILL range%2ld from%10.3e to%10.3eR in%4ld cell; ener res=%10.3e WIDTOT=%10.3e\n",
344 *ipnt,
345 rfield.anu[continuum.ifill0[*ipnt]] - rfield.widflx[continuum.ifill0[*ipnt]]/2.,
346 rfield.anu[continuum.ifill0[*ipnt]+nbin-1] + rfield.widflx[continuum.ifill0[*ipnt]+nbin-1]/2.,
347 nbin,
348 continuum.filres[*ipnt],
349 widtot );
350
351 fprintf( ioQQQ, " The requested range was%10.3e%10.3e The requested resolution was%10.3e\n",
352 fenlo, fenhi, resolv );
353 }
354
355 /* nrange is number of ranges */
356 *ipnt += 1;
357 continuum.nrange = MAX2(continuum.nrange,*ipnt);
358 return;
359}
360
361/*ChckFill perform sanity check confirming that the energy array has been properly filled */
362STATIC void ChckFill(void)
363{
364 bool lgFail;
365 long int i,
366 ipnt;
367 double energy;
368
369 DEBUG_ENTRY( "ChckFill()" );
370
371 ASSERT( rfield.anu[0] >= rfield.emm*0.99 );
372 ASSERT( rfield.anu[rfield.nupper-1] <= rfield.egamry*1.01 );
373
374 lgFail = false;
375 for( i=0; i < continuum.nrange; i++ )
376 {
377 /* test middle of energy bound */
378 energy = (continuum.filbnd[i] + continuum.filbnd[i+1])/2.;
379 ipnt = ipoint(energy);
380 if( energy < rfield.anu[ipnt-1] - rfield.widflx[ipnt-1]*0.5 )
381 {
382 fprintf( ioQQQ, " ChckFill middle test low fail\n" );
383 lgFail = true;
384 }
385
386 /* >>chng 02 jul 16, add second test - when "set resol 10" used,
387 * very large values of cell width, combined with fact that cells
388 * are log increasing, causes problem. */
389 else if( (energy > rfield.anu[ipnt-1] + rfield.widflx[ipnt-1]*0.5) &&
390 ( energy > rfield.anu[ipnt] - rfield.widflx[ipnt]*0.5 ) )
391 {
392 fprintf( ioQQQ, " ChckFill middle test high fail\n" );
393 lgFail = true;
394 }
395
396 /* test near low bound */
397 energy = continuum.filbnd[i]*0.99 + continuum.filbnd[i+1]*0.01;
398 ipnt = ipoint(energy);
399 if( energy < rfield.anu[ipnt-1] - rfield.widflx[ipnt-1]*0.5 )
400 {
401 fprintf( ioQQQ, " ChckFill low test low fail\n" );
402 lgFail = true;
403 }
404
405 else if( energy > rfield.anu[ipnt-1] + rfield.widflx[ipnt-1]* 0.5 )
406 {
407 fprintf( ioQQQ, " ChckFill low test high fail\n" );
408 lgFail = true;
409 }
410
411 /* test near high bound */
412 energy = continuum.filbnd[i]*0.01 + continuum.filbnd[i+1]*0.99;
413 ipnt = ipoint(energy);
414
415 if( energy < rfield.anu[ipnt-1] - rfield.widflx[ipnt-1]*0.5 )
416 {
417 fprintf( ioQQQ, " ChckFill high test low fail\n" );
418 lgFail = true;
419 }
420 /* >>chng 02 jul 16, add second test - when "set resol 10" used,
421 * very large values of cell width, combined with fact that cells
422 * are log increasing, causes problem. */
423 else if( (energy > rfield.anu[ipnt-1] + rfield.widflx[ipnt-1]*0.5) &&
424 ( energy > rfield.anu[ipnt] - rfield.widflx[ipnt]*0.5 ) )
425 {
426 fprintf( ioQQQ, " ChckFill high test high fail\n" );
427 lgFail = true;
428 }
429 }
430
431 if( lgFail )
432 {
434 }
435 return;
436}
437
438/* MALLOC arrays within rfield */
440{
441 long i;
442
443 DEBUG_ENTRY( "rfield_opac_malloc()" );
444
445 /* allocate one more than we use for the unit integration,
446 * will back up at end of routine */
447 ++rfield.nupper;
448
449 /* >>chng 03 feb 12, add fine mesh fine grid fine opacity array to keep track of line overlap */
453
454 /* frequency range in Rydberg needed for all resonance lines */
455 rfield.fine_ener_lo = rfield.emm;
456 rfield.fine_ener_hi = 1500.f;
457
458 /* set resolution of fine continuum mesh.
459 * rfield.fine_opac_velocity_width is width per cell, cm/s
460 * choose width so that most massive species (usually Fe) is well resolved
461 *
462 * rfield.fine_opac_nelem is the most massive (hence sharpest line)
463 * we will worry about. By default this is iron but can be changed
464 * with SET FINE CONTINUUM command
465 *
466 * TeLowestFineOpacity of 1e4 K is temperature were line width is
467 * evaluated. Tests were done using the stop temperature in its place
468 * Te below 1e4 K made fine opacity grid huge
469 * do not let temp get higher than 1e4 either - code run with stop temp 10 set
470 * stop temp of 1e10K and assert thrown at line 204 of cont_createpointers.c
471 * simply use 1e4 K as a characteristic temperature */
474 double TeLowestFineOpacity = 1e4;
475 rfield.fine_opac_velocity_width =
476 (realnum)sqrt(2.*BOLTZMANN/ATOMIC_MASS_UNIT*TeLowestFineOpacity/
477 dense.AtomicWeight[rfield.fine_opac_nelem] ) /
478 /* we want fine_opac_nresolv continuum elements across this line
479 * default is 1, changed with SET FINE CONTINUUM command */
480 rfield.fine_opac_nresolv;
481
482 /* we are at first zone so velocity shift is zero */
483 rfield.ipFineConVelShift = 0;
484
485 /* dimensionless resolution, dE/E, this is used in ipoint to get offset in find mesh */
486 rfield.fine_resol = rfield.fine_opac_velocity_width / SPEEDLIGHT;
487
488 /* the number of cells needed */
489 rfield.nfine_malloc = (long)(log10( rfield.fine_ener_hi / rfield.fine_ener_lo ) / log10( 1. + rfield.fine_resol ) );
490 if( rfield.nfine_malloc <= 0 )
492 rfield.nfine = rfield.nfine_malloc;
493
494 /* this is the fine opacity array to ghost the main low-resolution array */
495 rfield.fine_opac_zone = (realnum *)MALLOC(sizeof(realnum)*(unsigned)rfield.nfine_malloc );
496 memset(rfield.fine_opac_zone , 0 , (unsigned long)rfield.nfine_malloc*sizeof(realnum) );
497
498 /* this is the fine total optical array to ghost the main low-resolution array */
499 rfield.fine_opt_depth = (realnum *)MALLOC(sizeof(realnum)*(unsigned)rfield.nfine_malloc );
500 memset(rfield.fine_opt_depth , 0 , (unsigned long)rfield.nfine_malloc*sizeof(realnum) );
501
502 rfield.fine_anu = (realnum *)MALLOC(sizeof(realnum)*(unsigned)rfield.nfine_malloc );
503
504 /* now fill in energy array */
505 ASSERT( rfield.fine_ener_lo > 0. && rfield.fine_resol > 0 );
506 for( i=0;i<rfield.nfine_malloc; ++i )
507 {
508 rfield.fine_anu[i] = rfield.fine_ener_lo * (realnum)pow( (1.+rfield.fine_resol), (i+1.) );
509 }
510 /* done with fine array */
511
512 /* used to count number of lines per cell */
513 rfield.line_count = (long *)MALLOC(sizeof(long)*(unsigned)NCELL );
514 for( i=0; i<rfield.nupper; ++i)
515 {
516 rfield.line_count[i] = 0;
517 }
518 rfield.anu = (double*)MALLOC((size_t)(rfield.nupper*sizeof(double)) );
519 rfield.AnuOrg = (double*)MALLOC((size_t)(rfield.nupper*sizeof(double)) );
520 rfield.widflx = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
521 rfield.anulog = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
522 rfield.anusqr = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
523 rfield.anu2 = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
524 rfield.anu3 = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
525 rfield.flux_beam_time = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
526 rfield.flux_isotropic = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
527 rfield.flux_beam_const = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
528 rfield.flux_accum = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
529 rfield.ExtinguishFactor = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
530 rfield.convoc = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
531 rfield.OccNumbBremsCont = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
532 rfield.OccNumbIncidCont = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
533 rfield.OccNumbDiffCont = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
534 rfield.OccNumbContEmitOut = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
535 rfield.ConInterOut = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
536 rfield.SummedCon = (double*)MALLOC((size_t)(rfield.nupper*sizeof(double)) );
537 rfield.SummedDif = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
538 rfield.SummedDifSave = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
539 rfield.SummedOcc = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
540 rfield.ConOTS_local_photons = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
541 rfield.DiffuseEscape = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
542 rfield.TotDiff2Pht = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
543 rfield.ConOTS_local_OTS_rate = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
544 rfield.otslin = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
545 rfield.otscon = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
546 rfield.outlin_noplot = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
547 rfield.flux_beam_const_save = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
548 rfield.flux_time_beam_save = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
549 rfield.flux_isotropic_save = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
550 rfield.setCoarseTransCoefPtr((realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) ));
551 rfield.DiffuseLineEmission = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
552 rfield.ipnt_coarse_2_fine = (long int*)MALLOC((size_t)(rfield.nupper*sizeof(long int)) );
553
554 /* possibly save cumulative flux */
555 rfield.flux = (realnum**)MALLOC((size_t)(2*sizeof(realnum*)) );
556 rfield.ConEmitReflec = (realnum**)MALLOC((size_t)(2*sizeof(realnum*)) );
557 rfield.ConEmitOut = (realnum**)MALLOC((size_t)(2*sizeof(realnum*)) );
558 rfield.ConRefIncid = (realnum**)MALLOC((size_t)(2*sizeof(realnum*)) );
559 rfield.flux_total_incident = (realnum**)MALLOC((size_t)(2*sizeof(realnum*)) );
560 rfield.reflin = (realnum**)MALLOC((size_t)(2*sizeof(realnum*)) );
561 rfield.outlin = (realnum**)MALLOC((size_t)(2*sizeof(realnum*)) );
562
563 for( i=0; i<2; ++i )
564 {
565 rfield.flux[i] = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
566 rfield.ConEmitReflec[i] = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
567 rfield.ConEmitOut[i] = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
568 rfield.ConRefIncid[i] = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
569 rfield.flux_total_incident[i] = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
570 rfield.reflin[i] = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
571 rfield.outlin[i] = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
572 }
573 // the cumulative (time integral) emission
574 memset(rfield.flux[1] , 0 , (unsigned long)rfield.nupper*sizeof(realnum) );
575 memset(rfield.ConEmitReflec[1] , 0 , (unsigned long)rfield.nupper*sizeof(realnum) );
576 memset(rfield.ConEmitOut[1] , 0 , (unsigned long)rfield.nupper*sizeof(realnum) );
577 memset(rfield.ConRefIncid[1] , 0 , (unsigned long)rfield.nupper*sizeof(realnum) );
578 memset(rfield.flux_total_incident[1] , 0 , (unsigned long)rfield.nupper*sizeof(realnum) );
579 memset(rfield.reflin[1] , 0 , (unsigned long)rfield.nupper*sizeof(realnum) );
580 memset(rfield.outlin[1] , 0 , (unsigned long)rfield.nupper*sizeof(realnum) );
581
582 /* chng 02 may 16, by Ryan...added array for gaunt factors for ALL charges, malloc here. */
583 /* First index is EFFECTIVE CHARGE MINUS ONE! */
584 rfield.gff = (realnum**)MALLOC((size_t)((LIMELM+1)*sizeof(realnum*)) );
585 for( i = 1; i <= LIMELM; i++ )
586 {
587 rfield.gff[i] = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
588 }
589
590 rfield.csigh = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
591 rfield.csigc = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
592
593 rfield.comdn = (double*)MALLOC((size_t)(rfield.nupper*sizeof(double)) );
594 rfield.comup = (double*)MALLOC((size_t)(rfield.nupper*sizeof(double)) );
595 rfield.ContBoltz = (double*)MALLOC((size_t)(rfield.nupper*sizeof(double)) );
596
597 /*realnum rfield.otssav[NC_ELL][2];*/
598 rfield.otssav = (realnum**)MALLOC((size_t)(rfield.nupper*sizeof(realnum*)));
599 for( i=0; i<rfield.nupper; ++i)
600 {
601 rfield.otssav[i] = (realnum*)MALLOC(2*sizeof(realnum));
602 }
603
604
605 /* char rfield.chLineLabel[NLINES][5];*/
606 rfield.chLineLabel = (char**)MALLOC((size_t)(rfield.nupper*sizeof(char*)));
607 rfield.chContLabel = (char**)MALLOC((size_t)(rfield.nupper*sizeof(char*)));
608
609 /* now allocate all the labels for each of the above lines */
610 for( i=0; i<rfield.nupper; ++i)
611 {
612 rfield.chLineLabel[i] = (char*)MALLOC(5*sizeof(char));
613 rfield.chContLabel[i] = (char*)MALLOC(5*sizeof(char));
614 }
615
616 opac.TauAbsFace = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
617 memset( opac.TauAbsFace , 0 , rfield.nupper*sizeof(realnum) );
618
619 opac.TauScatFace = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
620 opac.E2TauAbsFace = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
621 opac.E2TauAbsTotal = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
622 opac.TauAbsTotal = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
623 opac.E2TauAbsOut = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
624 opac.ExpmTau = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
625 opac.tmn = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
626
627 opac.opacity_abs = (double*)MALLOC((size_t)(rfield.nupper*sizeof(double)) );
628 opac.opacity_abs_savzon1 = (double*)MALLOC((size_t)(rfield.nupper*sizeof(double)) );
629 opac.OldOpacSave = (double*)MALLOC((size_t)(rfield.nupper*sizeof(double)) );
630 opac.opacity_sct = (double*)MALLOC((size_t)(rfield.nupper*sizeof(double)) );
631 opac.albedo = (double*)MALLOC((size_t)(rfield.nupper*sizeof(double)) );
632 opac.opacity_sct_savzon1 = (double*)MALLOC((size_t)(rfield.nupper*sizeof(double)) );
633 opac.OpacStatic = (double*)MALLOC((size_t)(rfield.nupper*sizeof(double)) );
634 opac.FreeFreeOpacity = (double*)MALLOC((size_t)(rfield.nupper*sizeof(double)) );
635 opac.ExpZone = (double*)MALLOC((size_t)(rfield.nupper*sizeof(double)) );
636
637 opac.TauAbsGeo = (realnum**)MALLOC((size_t)(2*sizeof(realnum *)) );
638 opac.TauScatGeo = (realnum**)MALLOC((size_t)(2*sizeof(realnum *)) );
639 opac.TauTotalGeo = (realnum**)MALLOC((size_t)(2*sizeof(realnum *)) );
640
641 for( i=0; i<2; ++i)
642 {
643 opac.TauAbsGeo[i] = (realnum*)MALLOC(rfield.nupper*sizeof(realnum));
644 opac.TauScatGeo[i] = (realnum*)MALLOC(rfield.nupper*sizeof(realnum));
645 opac.TauTotalGeo[i] = (realnum*)MALLOC(rfield.nupper*sizeof(realnum));
646 }
647
648 /* fix allocate trick for one more than we use for the unit integration */
649 --rfield.nupper;
650
651 /* say that space exists */
652 lgRfieldMalloced = true;
653 return;
654}
655
656
657/* read the continuum definition from the file continuum_mesh.ini */
659{
660 FILE *ioDATA;
661 char chLine[INPUT_LINE_LENGTH];
662 long i;
663 bool lgEOL;
664 long i1 , i2 , i3;
665
666 DEBUG_ENTRY( "read_continuum_mesh()" );
667
668 if( trace.lgTrace )
669 fprintf( ioQQQ," read_continuum_mesh opening continuum_mesh.ini:");
670
671 ioDATA = open_data( "continuum_mesh.ini", "r" );
672
673 /* first line is a version number and does not count */
674 if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
675 {
676 fprintf( ioQQQ, " read_continuum_mesh could not read first line of continuum_mesh.ini.\n");
678 }
679 /* count how many lines are in the file, ignoring all lines
680 * starting with '#' */
681 continuum.nStoredBands = 0;
682 while( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) != NULL )
683 {
684 /* we want to count the lines that do not start with #
685 * since these contain data */
686 if( chLine[0] != '#')
687 ++continuum.nStoredBands;
688 }
689
690 /* we now have number of lines containing pairs of bounds,
691 * allocate space for the arrays we will need */
692 continuum.filbnd =
693 ((realnum *)MALLOC( (size_t)(continuum.nStoredBands+1)*sizeof(realnum )));
694 continuum.fildel =
695 ((realnum *)MALLOC( (size_t)(continuum.nStoredBands+1)*sizeof(realnum )));
696 continuum.filres =
697 ((realnum *)MALLOC( (size_t)(continuum.nStoredBands+1)*sizeof(realnum )));
698 continuum.ifill0 =
699 ((long *)MALLOC( (size_t)(continuum.nStoredBands+1)*sizeof(long )));
700 continuum.StoredEnergy =
701 ((double *)MALLOC( (size_t)(continuum.nStoredBands+1)*sizeof(double )));
702 continuum.StoredResolution =
703 ((double *)MALLOC( (size_t)(continuum.nStoredBands+1)*sizeof(double )));
704
705 /* now rewind the file so we can read it a second time*/
706 if( fseek( ioDATA , 0 , SEEK_SET ) != 0 )
707 {
708 fprintf( ioQQQ, " read_continuum_mesh could not rewind continuum_mesh.ini.\n");
710 }
711
712 /* check that magic number is ok */
713 if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
714 {
715 fprintf( ioQQQ, " read_continuum_mesh could not read first line of continuum_mesh.ini.\n");
717 }
718
719 i = 1;
720 /* continuum mesh magic number */
721 i1 = (long)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
722 i2 = (long)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
723 i3 = (long)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
724
725 bool lgResPower;
726
727 /* the following is the set of numbers that appear at the start of continuum_mesh.ini */
728 if( i1 == 1 && i2 == 9 && i3 == 29 )
729 // old version of the file (c08 and older), this has pairs: upper limit freq range, resolution
730 // this format is still supported to accomodate users with existing continuum_mesh.ini files.
731 lgResPower = false;
732 else if( i1 == 10 && i2 == 8 && i3 == 8 )
733 // new version of the file (c10 and newer), this has pairs: upper limit freq range, resolving power
734 // resolving power = 1./resolution
735 lgResPower = true;
736 else
737 {
738 fprintf( ioQQQ,
739 " read_continuum_mesh: the version of continuum_mesh.ini is not supported.\n" );
740 fprintf( ioQQQ,
741 " I found version number %li %li %li.\n" ,
742 i1 , i2 , i3 );
743 fprintf( ioQQQ, "Here is the line image:\n==%s==\n", chLine );
745 }
746
747 /* this starts at 1 not 0 since zero is reserved for the
748 * dummy line */
749 continuum.nStoredBands = 0;
750 while( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) != NULL )
751 {
752 /* only look at lines without '#' in first col */
753 if( chLine[0] != '#')
754 {
755 i = 1;
756 continuum.StoredEnergy[continuum.nStoredBands] = FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
757 continuum.StoredResolution[continuum.nStoredBands] = FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
758
759 // continuum energy could be 0 to indicate low or high energy bounds of code
760 // but none can be negative
761 if( continuum.StoredEnergy[continuum.nStoredBands]<0. ||
762 continuum.StoredResolution[continuum.nStoredBands]<=0. )
763 {
764 fprintf(ioQQQ, "DISASTER PROBLEM continuum_mesh.ini has a non-positive number.\n");
766 }
767
768 // convert resolving power (entered quantity) into resolution
769 if( lgResPower )
770 continuum.StoredResolution[continuum.nStoredBands] = 1./
771 continuum.StoredResolution[continuum.nStoredBands];
772
773 /* this is option to rescale resolution with set resolution command */
774 continuum.StoredResolution[continuum.nStoredBands] *= continuum.ResolutionScaleFactor;
775
776 ++continuum.nStoredBands;
777 }
778 }
779
780 fclose( ioDATA );
781
782 /* now verify continuum grid is ok - first are all values but the last positive? */
783 for( i=1; i<continuum.nStoredBands-1; ++i )
784 {
785 if( continuum.StoredEnergy[i-1]>=continuum.StoredEnergy[i] )
786 {
787 fprintf( ioQQQ,
788 " read_continuum_mesh: The continuum definition array energies must be in increasing order.\n" );
790 }
791 }
792 if( continuum.StoredEnergy[continuum.nStoredBands-1]!=0 )
793 {
794 fprintf( ioQQQ,
795 " read_continuum_mesh: The last continuum array energies must be zero.\n" );
797 }
798 return;
799}
800
801/*rfield_opac_zero zero out rfield arrays between certain limits */
803 /* index for first element in arrays to be set to zero */
804 long lo ,
805 /* array index for highest element to be set */
806 long ihi )
807{
808 long int i;
809
810 /* >>chng 01 aug 19, space not allocated yet,
811 * following code must also be present in contcreatemesh where
812 * space allocated for the first time */
813 if( lgRfieldMalloced )
814 {
815 unsigned long n=(unsigned long)(ihi-lo+1);
816 memset(&rfield.OccNumbDiffCont[lo] , 0 , n*sizeof(realnum) );
817 memset(&rfield.OccNumbContEmitOut[lo] , 0 , n*sizeof(realnum) );
818 memset(&rfield.ContBoltz[lo] , 0 , n*sizeof(double) );
819 /*>>chng 06 aug 15, this is now 2D array, saving diffuse continuum
820 * over all zones for use in exact RT */
821 /*memset(&rfield.ConEmitLocal[lo] , 0 , n*sizeof(realnum) );*/
822 memset(&rfield.ConEmitReflec[0][lo] , 0 , n*sizeof(realnum) );
823 memset(&rfield.ConEmitOut[0][lo] , 0 , n*sizeof(realnum) );
824 memset(&rfield.reflin[0][lo] , 0 , n*sizeof(realnum) );
825 memset(&rfield.ConRefIncid[0][lo] , 0 , n*sizeof(realnum) );
826 memset(&rfield.SummedCon[lo] , 0 , n*sizeof(double) );
827 memset(&rfield.OccNumbBremsCont[lo] , 0 , n*sizeof(realnum) );
828 memset(&rfield.convoc[lo] , 0 , n*sizeof(realnum) );
829 memset(&rfield.flux[0][lo] , 0 , n*sizeof(realnum) );
830 memset(&rfield.flux_total_incident[0][lo] , 0 , n*sizeof(realnum) );
831 memset(&rfield.flux_beam_const_save[lo] , 0 , n*sizeof(realnum) );
832 memset(&rfield.flux_time_beam_save[lo] , 0 , n*sizeof(realnum) );
833 memset(&rfield.flux_isotropic_save[lo] , 0 , n*sizeof(realnum) );
834 memset(&rfield.SummedOcc[lo] , 0 , n*sizeof(realnum) );
835 memset(&rfield.SummedDif[lo] , 0 , n*sizeof(realnum) );
836 memset(&rfield.flux_accum[lo] , 0 , n*sizeof(realnum) );
837 memset(&rfield.otslin[lo] , 0 , n*sizeof(realnum) );
838 memset(&rfield.otscon[lo] , 0 , n*sizeof(realnum) );
839 memset(&rfield.ConInterOut[lo] , 0 , n*sizeof(realnum) );
840 memset(&rfield.outlin[0][lo] , 0 , n*sizeof(realnum) );
841 memset(&rfield.outlin_noplot[lo] , 0 , n*sizeof(realnum) );
842 memset(&rfield.ConOTS_local_OTS_rate[lo], 0 , n*sizeof(realnum) );
843 memset(&rfield.ConOTS_local_photons[lo] , 0 , n*sizeof(realnum) );
844 memset(&opac.OldOpacSave[lo] , 0 , n*sizeof(double) );
845 memset(&opac.opacity_abs[lo] , 0 , n*sizeof(double) );
846 memset(&opac.opacity_sct[lo] , 0 , n*sizeof(double) );
847 memset(&opac.albedo[lo] , 0 , n*sizeof(double) );
848 memset(&opac.FreeFreeOpacity[lo] , 0 , n*sizeof(double) );
849
850 /* these are not defined on first iteration */
851 memset( &opac.E2TauAbsTotal[lo] , 0 , n*sizeof(realnum) );
852 memset( &opac.E2TauAbsOut[lo] , 0 , n*sizeof(realnum) );
853 memset( &opac.TauAbsTotal[lo] , 0 , n*sizeof(realnum) );
854
855 for( i=lo; i <= ihi; i++ )
856 {
857 opac.TauTotalGeo[0][i] = opac.taumin;
858 opac.TauAbsGeo[0][i] = opac.taumin;
859 opac.TauScatGeo[0][i] = opac.taumin;
860 opac.tmn[i] = 1.;
861 opac.ExpZone[i] = 1.;
862 opac.E2TauAbsFace[i] = 1.;
863 opac.ExpmTau[i] = 1.;
864 opac.OpacStatic[i] = 1.;
865 }
866 /* also zero out fine opacity fine grid fine mesh array */
867 memset(rfield.fine_opac_zone , 0 , (unsigned long)rfield.nfine_malloc*sizeof(realnum) );
868 /* also zero out fine opacity array */
869 memset(rfield.fine_opt_depth , 0 , (unsigned long)rfield.nfine_malloc*sizeof(realnum) );
870 }
871 return;
872}
FILE * ioQQQ
Definition cddefines.cpp:7
#define ASSERT(exp)
Definition cddefines.h:578
#define STATIC
Definition cddefines.h:97
#define MALLOC(exp)
Definition cddefines.h:501
const int LIMELM
Definition cddefines.h:258
const int INPUT_LINE_LENGTH
Definition cddefines.h:254
#define EXIT_FAILURE
Definition cddefines.h:140
double FFmtRead(const char *chCard, long int *ipnt, long int last, bool *lgEOL)
Definition service.cpp:381
#define cdEXIT(FAIL)
Definition cddefines.h:434
float realnum
Definition cddefines.h:103
#define MAX2
Definition cddefines.h:782
char * read_whole_line(char *chLine, int nChar, FILE *ioIN)
Definition service.cpp:70
NORETURN void TotalInsanity(void)
Definition service.cpp:886
#define DEBUG_ENTRY(funcname)
Definition cddefines.h:684
void ShowMe(void)
Definition service.cpp:181
bool lgRfieldMalloced
Definition cdinit.cpp:98
STATIC void rfield_opac_malloc(void)
STATIC void fill(double fenlo, double fenhi, double resolv, long int *n0, long int *ipnt, bool lgCount)
void rfield_opac_zero(long lo, long ihi)
void ContCreateMesh(void)
STATIC void read_continuum_mesh(void)
STATIC void ChckFill(void)
long ipoint(double energy_ryd)
t_continuum continuum
Definition continuum.cpp:5
FILE * open_data(const char *fname, const char *mode, access_scheme scheme)
Definition cpu.cpp:625
t_dense dense
Definition dense.cpp:24
t_geometry geometry
Definition geometry.cpp:5
t_iterations iterations
Definition iterations.cpp:5
t_opac opac
Definition opacity.cpp:5
UNUSED const double SPEEDLIGHT
Definition physconst.h:100
UNUSED const double BOLTZMANN
Definition physconst.h:97
UNUSED const double ATOMIC_MASS_UNIT
Definition physconst.h:88
t_rfield rfield
Definition rfield.cpp:8
const int NCELL
Definition rfield.h:21
t_trace trace
Definition trace.cpp:5