cloudy  trunk
 All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
cont_createmesh.cpp
Go to the documentation of this file.
1 /* This file is part of Cloudy and is copyright (C)1978-2008 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 "doppvel.h"
14 #include "dense.h"
15 #include "trace.h"
16 #include "opacity.h"
17 #include "ipoint.h"
18 #include "geometry.h"
19 #include "continuum.h"
20 
21 /* read the continuum definition from the file continuum_mesh.ini */
22 STATIC void read_continuum_mesh( void );
23 
24 /*fill define the continuum energy grid over a specified range */
25 STATIC void fill(double fenlo,
26  double fenhi,
27  double resolv,
28  long int *n0,
29  long int *ipnt,
30  /* this says only count, do not fill */
31  bool lgCount );
32 
33 /*rfield_opac_malloc MALLOC space for opacity arrays */
34 STATIC void rfield_opac_malloc(void);
35 
36 /*ChckFill perform sanity check confirming that the energy array has been properly filled */
37 STATIC void ChckFill(void);
38 
39 void ContCreateMesh(void)
40 {
41  long int
42  i,
43  ipnt,
44  n0;
45 
46  /* flag to say whether pointers have ever been evaluated */
47  static bool lgPntEval = false;
48 
49  DEBUG_ENTRY( "ContCreateMesh()" );
50 
51  /* lgPntEval is local static variable defined false when defined.
52  * it is set true below, so that pointers only created one time in the
53  * history of this coreload. */
54  if( lgPntEval )
55  {
56  if( trace.lgTrace )
57  {
58  fprintf( ioQQQ, " ContCreateMesh called, not evaluating.\n" );
59  }
60  /* now save current form of energy array */
61  for( i=0; i < rfield.nupper; i++ )
62  {
63  rfield.anu[i] = rfield.AnuOrg[i];
64  rfield.anu2[i] = rfield.anu[i]*rfield.anu[i];
65  }
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;
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] ,
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.");
128  cdEXIT(EXIT_FAILURE);
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 */
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 */
143  for( i=1; i < iterations.iter_malloc; i++ )
144  {
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  for( i=0;i<(geometry.nend_max+1); ++i )
150  {
151  rfield.ConEmitLocal[i] = (realnum*)MALLOC( (size_t)rfield.nupper*sizeof(realnum) );
152  }
153 
154  /* ================================================================ */
155  n0 = 1;
156  ipnt = 0;
157 
158  /* this is number of ranges that will be introduced below*/
159  continuum.nrange = 0;
160 
161  /* the default array values are set in continuum_mesh.ini */
163  &n0,&ipnt,false);
164  for(i=1; i<continuum.nStoredBands; ++i )
165  {
166  fill(continuum.StoredEnergy[i-1] ,
169  &n0,&ipnt,false);
170  }
171 
172  /* ================================================================ */
173 
174  /* fill in the false highest cell used for unit verification */
178 
179  /* there must be a cell above nflux for us to pass unity through the vol integrator
180  * as a sanity check. assert that this is true so we will crash if ever changed
181  ASSERT( rfield.nupper +1 <= rfield.nupper );*/
182 
183  /* this is done here when the space is first allocated,
184  * then done on every subsequent initialization in zero.c */
186 
187  /* this is a sanity check for results produced above by fill */
188  ChckFill();
189 
190  /* now fix widflx array so that it is correct */
191  for( i=1; i<rfield.nupper-1; ++i )
192  {
193  rfield.widflx[i] = ((rfield.anu[i+1] - rfield.anu[i]) + (rfield.anu[i] -
194  rfield.anu[i-1]))/2.f;
195  }
196 
197  ipnt = 0;
198  /* now save current form of array, and define some quantities related to it */
199  for( i=0; i < rfield.nupper; i++ )
200  {
201  double alf , bet;
202 
203  rfield.AnuOrg[i] = rfield.anu[i];
204  rfield.anusqr[i] = (realnum)sqrt(rfield.AnuOrg[i]);
205  /* following are Compton exchange factors from Tarter */
206  /* this code also appears in highen, but coef needed before that routine called. */
207  alf = 1./(1. + rfield.anu[i]*(1.1792e-4 + 7.084e-10*rfield.anu[i]));
208  bet = 1. - alf*rfield.anu[i]*(1.1792e-4 + 2.*7.084e-10*rfield.anu[i])/4.;
209  rfield.csigh[i] = (realnum)(alf*rfield.anu[i]*rfield.anu[i]*3.858e-25);
210  rfield.csigc[i] = (realnum)(alf*bet*rfield.anu[i]*3.858e-25);
211  rfield.anu2[i] = rfield.anu[i]*rfield.anu[i];
212 
213  /* >>chng 05 feb 28, add transmission and mapping coef */
214  /* map these coarse continua into fine continuum grid */
216  {
217  /* 0 (false) says not defined */
218  rfield.ipnt_coarse_2_fine[i] = 0;
219  }
220  else
221  {
222  if( ipnt==0 )
223  {
224  /* this is the first one that maps onto the fine array */
225  rfield.ipnt_coarse_2_fine[i] = 0;
226  ipnt = 1;
227  }
228  else
229  {
230  /* find first fine frequency that is greater than this coarse value */
231  while (ipnt < rfield.nfine_malloc && rfield.fine_anu[ipnt]<rfield.anu[i] )
232  {
233  ++ipnt;
234  }
235  rfield.ipnt_coarse_2_fine[i] = ipnt;
236  }
237  }
238  /*fprintf(ioQQQ," coarse %li nu= %.3e points to fine %li nu=%.3e\n",
239  i, rfield.anu[i] , rfield.ipnt_coarse_2_fine[i] , rfield.fine_anu[rfield.ipnt_coarse_2_fine[i]] );*/
240  rfield.trans_coef_zone[i] = 1.f;
241  rfield.trans_coef_total[i] = 1.f;
242  }
243  return;
244 }
245 
246 /*fill define the continuum energy grid over a specified range, called by ContCreateMesh */
248  /* lower bounds to this energy range */
249  double fenlo,
250  /* upper bounds to this continuum range */
251  double fenhi,
252  /* relative energy resolution */
253  double resolv,
254  /* starting index within frequency grid */
255  long int *n0,
256  /* which energy band this is */
257  long int *ipnt,
258  /* this says only count, do not fill */
259  bool lgCount )
260 {
261  long int i,
262  nbin;
263  realnum widtot;
264  double aaa , bbb;
265 
266  DEBUG_ENTRY( "fill()" );
267 
268  ASSERT( fenlo>0. && fenhi>0. && resolv>0. );
269 
270  /* this is the number of cells needed to fill the array with numbers at the requested resolution */
271  nbin = (long int)(log(10.)*log10(fenhi/fenlo)/resolv + 1);
272 
273  if( lgCount )
274  {
275  /* true means only count number of cells, don't do anything */
276  *n0 += nbin;
277  return;
278  }
279 
280  if( *ipnt > 0 && fabs(1.-fenlo/continuum.filbnd[*ipnt]) > 1e-4 )
281  {
282  fprintf( ioQQQ, " FILL improper bounds.\n" );
283  fprintf( ioQQQ, " ipnt=%3ld fenlo=%11.4e filbnd(ipnt)=%11.4e\n",
284  *ipnt, fenlo, continuum.filbnd[*ipnt] );
285  cdEXIT(EXIT_FAILURE);
286  }
287 
288  ASSERT( *ipnt < continuum.nStoredBands );
289 
290  continuum.ifill0[*ipnt] = *n0 - 1;
291  continuum.filbnd[*ipnt] = (realnum)fenlo;
292  continuum.filbnd[*ipnt+1] = (realnum)fenhi;
293 
294  /* this is the number of cells needed to fill the array with numbers
295  nbin = (long int)(log(10.)*log10(fenhi/fenlo)/resolv + 1);*/
296  continuum.fildel[*ipnt] = (realnum)(log10(fenhi/fenlo)/nbin);
297 
298  if( continuum.fildel[*ipnt] < 0.01 )
299  {
300  continuum.filres[*ipnt] = (realnum)(log(10.)*continuum.fildel[*ipnt]);
301  }
302  else
303  {
304  continuum.filres[*ipnt] = (realnum)((pow(10.,2.*continuum.fildel[*ipnt]) - 1.)/2./
305  pow((realnum)10.f,continuum.fildel[*ipnt]));
306  }
307 
308  if( (*n0 + nbin-2) > rfield.nupper )
309  {
310  fprintf( ioQQQ, " Fill would need %ld cells to get to an energy of %.3e\n",
311  *n0 + nbin, fenhi );
312  fprintf( ioQQQ, " This is a major logical error in fill.\n");
313  ShowMe();
314  cdEXIT(EXIT_FAILURE);
315  }
316 
317  widtot = 0.;
318  for( i=0; i < nbin; i++ )
319  {
320  bbb = continuum.fildel[*ipnt]*((realnum)(i) + 0.5);
321  aaa = pow( 10. , bbb );
322 
323  rfield.anu[i+continuum.ifill0[*ipnt]] = (realnum)(fenlo*aaa);
324 
325  rfield.widflx[i+continuum.ifill0[*ipnt]] = rfield.anu[i+continuum.ifill0[*ipnt]]*
326  continuum.filres[*ipnt];
327 
328  widtot += rfield.widflx[i+continuum.ifill0[*ipnt]];
329  }
330 
331  *n0 += nbin;
332  if( trace.lgTrace && (trace.lgConBug || trace.lgPtrace) )
333  {
334  fprintf( ioQQQ,
335  " FILL range%2ld from%10.3e to%10.3eR in%4ld cell; ener res=%10.3e WIDTOT=%10.3e\n",
336  *ipnt,
337  rfield.anu[continuum.ifill0[*ipnt]] - rfield.widflx[continuum.ifill0[*ipnt]]/2.,
338  rfield.anu[continuum.ifill0[*ipnt]+nbin-1] + rfield.widflx[continuum.ifill0[*ipnt]+nbin-1]/2.,
339  nbin,
340  continuum.filres[*ipnt],
341  widtot );
342 
343  fprintf( ioQQQ, " The requested range was%10.3e%10.3e The requested resolution was%10.3e\n",
344  fenlo, fenhi, resolv );
345  }
346 
347  /* nrange is number of ranges */
348  *ipnt += 1;
350  return;
351 }
352 
353 /*ChckFill perform sanity check confirming that the energy array has been properly filled */
354 STATIC void ChckFill(void)
355 {
356  bool lgFail;
357  long int i,
358  ipnt;
359  double energy;
360 
361  DEBUG_ENTRY( "ChckFill()" );
362 
363  ASSERT( rfield.anu[0] >= rfield.emm*0.99 );
364  ASSERT( rfield.anu[rfield.nupper-1] <= rfield.egamry*1.01 );
365 
366  lgFail = false;
367  for( i=0; i < continuum.nrange; i++ )
368  {
369  /* test middle of energy bound */
370  energy = (continuum.filbnd[i] + continuum.filbnd[i+1])/2.;
371  ipnt = ipoint(energy);
372  if( energy < rfield.anu[ipnt-1] - rfield.widflx[ipnt-1]*0.5 )
373  {
374  fprintf( ioQQQ, " ChckFill middle test low fail\n" );
375  lgFail = true;
376  }
377 
378  /* >>chng 02 jul 16, add second test - when "set resol 10" used,
379  * very large values of cell width, combined with fact that cells
380  * are log increasing, causes problem. */
381  else if( (energy > rfield.anu[ipnt-1] + rfield.widflx[ipnt-1]*0.5) &&
382  ( energy > rfield.anu[ipnt] - rfield.widflx[ipnt]*0.5 ) )
383  {
384  fprintf( ioQQQ, " ChckFill middle test high fail\n" );
385  lgFail = true;
386  }
387 
388  /* test near low bound */
389  energy = continuum.filbnd[i]*0.99 + continuum.filbnd[i+1]*0.01;
390  ipnt = ipoint(energy);
391  if( energy < rfield.anu[ipnt-1] - rfield.widflx[ipnt-1]*0.5 )
392  {
393  fprintf( ioQQQ, " ChckFill low test low fail\n" );
394  lgFail = true;
395  }
396 
397  else if( energy > rfield.anu[ipnt-1] + rfield.widflx[ipnt-1]* 0.5 )
398  {
399  fprintf( ioQQQ, " ChckFill low test high fail\n" );
400  lgFail = true;
401  }
402 
403  /* test near high bound */
404  energy = continuum.filbnd[i]*0.01 + continuum.filbnd[i+1]*0.99;
405  ipnt = ipoint(energy);
406 
407  if( energy < rfield.anu[ipnt-1] - rfield.widflx[ipnt-1]*0.5 )
408  {
409  fprintf( ioQQQ, " ChckFill high test low fail\n" );
410  lgFail = true;
411  }
412  /* >>chng 02 jul 16, add second test - when "set resol 10" used,
413  * very large values of cell width, combined with fact that cells
414  * are log increasing, causes problem. */
415  else if( (energy > rfield.anu[ipnt-1] + rfield.widflx[ipnt-1]*0.5) &&
416  ( energy > rfield.anu[ipnt] - rfield.widflx[ipnt]*0.5 ) )
417  {
418  fprintf( ioQQQ, " ChckFill high test high fail\n" );
419  lgFail = true;
420  }
421  }
422 
423  if( lgFail )
424  {
425  cdEXIT(EXIT_FAILURE);
426  }
427  return;
428 }
429 
430 /* MALLOC arrays within rfield */
432 {
433  long i;
434 
435  DEBUG_ENTRY( "rfield_opac_malloc()" );
436 
437  /* allocate one more than we use for the unit integration,
438  * will back up at end of routine */
439  ++rfield.nupper;
440 
441  /* >>chng 03 feb 12, add fine mesh fine grid fine opacity array to keep track of line overlap */
446  /* frequency range in Rydberg needed for all resonance lines */
447  rfield.fine_ener_lo = 0.05f;
448  rfield.fine_ener_hi = 1500.f;
449 
450  /* set resolution of fine continuum mesh.
451  * rfield.fine_opac_velocity_width is width per cell, cm/s
452  * choose width so that most massive species (usually Fe) is well resolved
453  *
454  * rfield.fine_opac_nelem is the most massive (hence sharpest line)
455  * we will worry about. By default this is irion but can be changed
456  * with SET FINE CONTINUUM command
457  *
458  * TeLowestFineOpacity of 1e4 K is temperature were line width is
459  * evaluated. Tests were done using the stop temperature in its place
460  * Te below 1e4 K made fine opacity grid huge
461  * do not let temp get higher than 1e4 either - code run with stop temp 10 set
462  * stop temp of 1e10K and assert thrown at line 204 of cont_createpointers.c
463  * simply use 1e4 K as a characteristic temperature */
466  double TeLowestFineOpacity = 1e4;
468  (realnum)sqrt(2.*BOLTZMANN/ATOMIC_MASS_UNIT*TeLowestFineOpacity/
470  /* we want fine_opac_nresolv continuum elements across this line
471  * default is 1, changed with SET FINE CONTINUUM command */
473 
474  /* we are at first zone so velocity shift is zero */
476 
477  /* dimensionless resolution, dE/E, this is used in ipoint to get offset in find mesh */
479 
480  /* the number of cells needed */
481  rfield.nfine_malloc = (long)(log10( rfield.fine_ener_hi / rfield.fine_ener_lo ) / log10( 1. + rfield.fine_resol ) );
482  if( rfield.nfine_malloc <= 0 )
483  TotalInsanity();
485 
486  /* this is the fine opacity array to ghost the main low-resolution array */
487  rfield.fine_opac_zone = (realnum *)MALLOC(sizeof(realnum)*(unsigned)rfield.nfine_malloc );
488  memset(rfield.fine_opac_zone , 0 , (unsigned long)rfield.nfine_malloc*sizeof(realnum) );
489 
490  /* this is the fine total optical array to ghost the main low-resolution array */
491  rfield.fine_opt_depth = (realnum *)MALLOC(sizeof(realnum)*(unsigned)rfield.nfine_malloc );
492  memset(rfield.fine_opt_depth , 0 , (unsigned long)rfield.nfine_malloc*sizeof(realnum) );
493 
494  rfield.fine_anu = (realnum *)MALLOC(sizeof(realnum)*(unsigned)rfield.nfine_malloc );
495 
496  /* now fill in energy array */
497  ASSERT( rfield.fine_ener_lo > 0. && rfield.fine_resol > 0 );
498  for( i=0;i<rfield.nfine_malloc; ++i )
499  {
500  rfield.fine_anu[i] = rfield.fine_ener_lo * (realnum)pow( (1.+rfield.fine_resol), (i+1.) );
501  }
502  /* done with fine array */
503 
504  /* used to count number of lines per cell */
505  rfield.line_count = (long *)MALLOC(sizeof(long)*(unsigned)NCELL );
506  for( i=0; i<rfield.nupper; ++i)
507  {
508  rfield.line_count[i] = 0;
509  }
510  rfield.anu = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
511  rfield.AnuOrg = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
512  rfield.widflx = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
513  rfield.anulog = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
514  rfield.anusqr = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
515  rfield.anu2 = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
516  rfield.anu3 = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
517  rfield.flux_beam_time = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
518  rfield.flux_isotropic = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
519  rfield.flux_beam_const = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
520  rfield.flux = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
521  rfield.flux_accum = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
522  rfield.convoc = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
523  rfield.OccNumbBremsCont = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
524  rfield.OccNumbIncidCont = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
525  rfield.OccNumbDiffCont = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
526  rfield.OccNumbContEmitOut = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
527  rfield.ConEmitReflec = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
528  rfield.ConEmitOut = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
529  rfield.ConInterOut = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
530  rfield.ConRefIncid = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
531  rfield.SummedCon = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
532  rfield.SummedDif = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
533  rfield.SummedDifSave = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
534  rfield.SummedOcc = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
536  rfield.TotDiff2Pht = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
538  rfield.otslin = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
539  rfield.otscon = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
540  rfield.otslinNew = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
541  rfield.otsconNew = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
542  rfield.outlin = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
543  rfield.outlin_noplot = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
544  rfield.reflin = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
545  rfield.flux_total_incident = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
547  rfield.flux_time_beam_save = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
548  rfield.flux_isotropic_save = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
549  rfield.trans_coef_zone = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
550  rfield.trans_coef_total = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
551  rfield.ipnt_coarse_2_fine = (long int*)MALLOC((size_t)(rfield.nupper*sizeof(long int)) );
552 
553  /* chng 02 may 16, by Ryan...added array for gaunt factors for ALL charges, malloc here. */
554  /* First index is EFFECTIVE CHARGE MINUS ONE! */
555  rfield.gff = (realnum**)MALLOC((size_t)((LIMELM+1)*sizeof(realnum*)) );
556  for( i = 1; i <= LIMELM; i++ )
557  {
558  rfield.gff[i] = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
559  }
560 
561  rfield.csigh = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
562  rfield.csigc = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
563  rfield.ConTabRead = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
564 
565  rfield.comdn = (double*)MALLOC((size_t)(rfield.nupper*sizeof(double)) );
566  rfield.comup = (double*)MALLOC((size_t)(rfield.nupper*sizeof(double)) );
567  rfield.ContBoltz = (double*)MALLOC((size_t)(rfield.nupper*sizeof(double)) );
568 
569  /*realnum rfield.otssav[NC_ELL][2];*/
570  rfield.otssav = (realnum**)MALLOC((size_t)(rfield.nupper*sizeof(realnum*)));
571  for( i=0; i<rfield.nupper; ++i)
572  {
573  rfield.otssav[i] = (realnum*)MALLOC(2*sizeof(realnum));
574  }
575 
576 
577  /* char rfield.chLineLabel[NLINES][5];*/
578  rfield.chLineLabel = (char**)MALLOC((size_t)(rfield.nupper*sizeof(char*)));
579  rfield.chContLabel = (char**)MALLOC((size_t)(rfield.nupper*sizeof(char*)));
580 
581  /* now allocate all the labels for each of the above lines */
582  for( i=0; i<rfield.nupper; ++i)
583  {
584  rfield.chLineLabel[i] = (char*)MALLOC(5*sizeof(char));
585  rfield.chContLabel[i] = (char*)MALLOC(5*sizeof(char));
586  }
587 
588  opac.TauAbsFace = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
589  memset( opac.TauAbsFace , 0 , rfield.nupper*sizeof(realnum) );
590 
591  opac.TauScatFace = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
592  opac.E2TauAbsFace = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
593  opac.E2TauAbsTotal = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
594  opac.TauAbsTotal = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
595  opac.E2TauAbsOut = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
596  opac.ExpmTau = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
597  opac.tmn = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
598 
599  opac.opacity_abs = (double*)MALLOC((size_t)(rfield.nupper*sizeof(double)) );
600  opac.opacity_abs_savzon1 = (double*)MALLOC((size_t)(rfield.nupper*sizeof(double)) );
601  opac.OldOpacSave = (double*)MALLOC((size_t)(rfield.nupper*sizeof(double)) );
602  opac.opacity_sct = (double*)MALLOC((size_t)(rfield.nupper*sizeof(double)) );
603  opac.albedo = (double*)MALLOC((size_t)(rfield.nupper*sizeof(double)) );
604  opac.opacity_sct_savzon1 = (double*)MALLOC((size_t)(rfield.nupper*sizeof(double)) );
605  opac.OpacStatic = (double*)MALLOC((size_t)(rfield.nupper*sizeof(double)) );
606  opac.FreeFreeOpacity = (double*)MALLOC((size_t)(rfield.nupper*sizeof(double)) );
607  opac.ExpZone = (double*)MALLOC((size_t)(rfield.nupper*sizeof(double)) );
608 
609  opac.TauAbsGeo = (realnum**)MALLOC((size_t)(2*sizeof(realnum *)) );
610  opac.TauScatGeo = (realnum**)MALLOC((size_t)(2*sizeof(realnum *)) );
611  opac.TauTotalGeo = (realnum**)MALLOC((size_t)(2*sizeof(realnum *)) );
612 
613  for( i=0; i<2; ++i)
614  {
615  opac.TauAbsGeo[i] = (realnum*)MALLOC(rfield.nupper*sizeof(realnum));
618  }
619 
620  /* fix allocate trick for one more than we use for the unit integration */
621  --rfield.nupper;
622 
623  /* say that space exists */
624  lgRfieldMalloced = true;
625  return;
626 }
627 
628 
629 /* read the continuum definition from the file continuum_mesh.ini */
631 {
632  FILE *ioDATA;
633  char chLine[INPUT_LINE_LENGTH];
634  long i;
635  bool lgEOL;
636  long i1 , i2 , i3;
637 
638  DEBUG_ENTRY( "read_continuum_mesh()" );
639 
640  if( trace.lgTrace )
641  fprintf( ioQQQ," read_continuum_mesh opening continuum_mesh.ini:");
642 
643  ioDATA = open_data( "continuum_mesh.ini", "r" );
644 
645  /* first line is a version number and does not count */
646  if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
647  {
648  fprintf( ioQQQ, " read_continuum_mesh could not read first line of continuum_mesh.ini.\n");
649  cdEXIT(EXIT_FAILURE);
650  }
651  /* count how many lines are in the file, ignoring all lines
652  * starting with '#' */
654  while( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) != NULL )
655  {
656  /* we want to count the lines that do not start with #
657  * since these contain data */
658  if( chLine[0] != '#')
660  }
661 
662  /* we now have number of lines containing pairs of bounds,
663  * allocate space for the arrays we will need */
664  continuum.filbnd =
665  ((realnum *)MALLOC( (size_t)(continuum.nStoredBands+1)*sizeof(realnum )));
666  continuum.fildel =
667  ((realnum *)MALLOC( (size_t)(continuum.nStoredBands+1)*sizeof(realnum )));
668  continuum.filres =
669  ((realnum *)MALLOC( (size_t)(continuum.nStoredBands+1)*sizeof(realnum )));
670  continuum.ifill0 =
671  ((long *)MALLOC( (size_t)(continuum.nStoredBands+1)*sizeof(long )));
673  ((double *)MALLOC( (size_t)(continuum.nStoredBands+1)*sizeof(double )));
675  ((double *)MALLOC( (size_t)(continuum.nStoredBands+1)*sizeof(double )));
676 
677  /* now rewind the file so we can read it a second time*/
678  if( fseek( ioDATA , 0 , SEEK_SET ) != 0 )
679  {
680  fprintf( ioQQQ, " read_continuum_mesh could not rewind continuum_mesh.ini.\n");
681  cdEXIT(EXIT_FAILURE);
682  }
683 
684  /* check that magic number is ok */
685  if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
686  {
687  fprintf( ioQQQ, " read_continuum_mesh could not read first line of continuum_mesh.ini.\n");
688  cdEXIT(EXIT_FAILURE);
689  }
690 
691  i = 1;
692  /* level 1 magic number */
693  i1 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
694  i2 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
695  i3 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
696 
697  /* the following is the set of numbers that appear at the start of continuum_mesh.ini 01 08 10 */
698  if( ( i1 != 1 ) || ( i2 != 9 ) || ( i3 != 29 ) )
699  {
700  fprintf( ioQQQ,
701  " read_continuum_mesh: the version of continuum_mesh.ini is not the current version.\n" );
702  fprintf( ioQQQ,
703  " I expected to find the number 01 09 29 and got %li %li %li instead.\n" ,
704  i1 , i2 , i3 );
705  fprintf( ioQQQ, "Here is the line image:\n==%s==\n", chLine );
706  cdEXIT(EXIT_FAILURE);
707  }
708 
709  /* this starts at 1 not 0 since zero is reserved for the
710  * dummy line */
712  while( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) != NULL )
713  {
714  /* only look at lines without '#' in first col */
715  if( chLine[0] != '#')
716  {
717  i = 1;
720 
721  /* option to enter numbers as logs if less than zero */
725 
729 
730  /* this is option to rescale resolution with set resolution command */
732 
734  {
735  fprintf( ioQQQ,
736  " read_continuum_mesh: A continuum resolution was zero - this is not allowed.\n" );
737  cdEXIT(EXIT_FAILURE);
738  }
739 
741  }
742  }
743 
744  fclose( ioDATA );
745 
746  /* now verify continuum grid is ok - first are all values but the last positive? */
747  for( i=1; i<continuum.nStoredBands-1; ++i )
748  {
750  {
751  fprintf( ioQQQ,
752  " read_continuum_mesh: The continuum definition array energies must be in increasing order.\n" );
753  cdEXIT(EXIT_FAILURE);
754  }
755  }
757  {
758  fprintf( ioQQQ,
759  " read_continuum_mesh: The last continuum array energies must be zero.\n" );
760  cdEXIT(EXIT_FAILURE);
761  }
762  return;
763 }

Generated for cloudy by doxygen 1.8.1.2