cloudy  trunk
 All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
parse_drive.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 /*ParseDrive parse the drive command - drive calls to various subs */
4 /*DrvCaseBHS allow user to query hydrogen A's, asks for up, low level, gives A, drive hyas */
5 /*DrvHyas allow user to query hydrogen A's, asks for up, low level, gives A, drive hyas */
6 /*dgaunt drive gaunt factor routines by letting user query values */
7 #include "cddefines.h"
8 #include "trace.h"
9 #include "hydro_bauman.h"
10 #include "hydroeinsta.h"
11 #include "phycon.h"
12 #include "atmdat.h"
13 #include "taulines.h"
14 #include "thermal.h"
15 #include "abund.h"
16 #include "rt.h"
17 #include "continuum.h"
18 #include "parse.h"
19 #include "physconst.h"
20 
21 /*dgaunt drive gaunt factor routines by letting user query values */
22 STATIC void dgaunt(void);
23 
24 /*DrvHyas allow user to query hydrogen A's, asks for up, low level, gives A, drive hyas */
25 STATIC void DrvHyas(void);
26 
27 /* drive escape probability routines */
28 STATIC void DrvEscP( void );
29 
30 /*DrvCaseBHS allow user to query hydrogen A's, asks for up, low level, gives A, drive hyas */
31 STATIC void DrvCaseBHS(void);
32 
33 void ParseDrive(char *chCard )
34 {
35  bool lgEOL;
36  long int n,
37  i;
38  double fac,
39  zed;
40 
41  DEBUG_ENTRY( "ParseDrive()" );
42 
43  /* NB evolve all following names to style DrvSomething */
44 
45  /* option to drive cloudy, which one? */
46  if( nMatch("FFMT",chCard) )
47  {
48  /* free format parser */
49  char chInput[INPUT_LINE_LENGTH];
50  fprintf( ioQQQ, " FFmtRead ParseDrive entered. Enter number.\n" );
51  lgEOL = false;
52  while( !lgEOL )
53  {
54  if( read_whole_line( chInput , (int)sizeof(chInput) , ioStdin ) == NULL )
55  {
56  fprintf( ioQQQ, " ParseDrive.dat error getting magic number\n" );
57  cdEXIT(EXIT_FAILURE);
58  }
59  i = 1;
60  fac = FFmtRead(chInput,&i,INPUT_LINE_LENGTH,&lgEOL);
61  if( lgEOL )
62  {
63  fprintf( ioQQQ, " FFmtRead hit the EOL with no value, return=%10.2e\n",
64  fac );
65  break;
66  }
67  else if( fac == 0. )
68  {
69  break;
70  }
71  else
72  {
73  fprintf( ioQQQ, " FFmtRead returned with value%11.4e\n",
74  fac );
75  }
76  fprintf( ioQQQ, " Enter 0 to stop, or another value.\n" );
77  }
78  fprintf( ioQQQ, " FFmtRead ParseDrive exits.\n" );
79  }
80 
81  else if( nMatch("CASE",chCard) )
82  {
83  /* option to interpolate on Hummer and Storey case b hydrogenic emission routines */
84  DrvCaseBHS( );
85  }
86 
87  else if( nMatch("CDLI",chCard) )
88  {
89  /* drive cdLine to check that it finds all the right lines, routine is in lines.c */
90  trace.lgDrv_cdLine = true;
91  }
92 
93  else if( nMatch(" E1 ",chCard) )
94  {
95  /* option to drive exponential integral routines */
96  for( i=0; i<30; ++i )
97  {
98  double tau = pow( 10. , ((double)i/4. - 5.) );
99  fprintf(ioQQQ,"tau\t%.3e\t exp-tau\t%.5e\t e1 tau\t%.5e \t ee2 tau\t%.5e \n",
100  tau , sexp(tau) , ee1(tau) , e2(tau ) );
101  }
102  cdEXIT(EXIT_SUCCESS);
103  }
104 
105  else if( nMatch("ESCA",chCard) )
106  {
107  /* option to drive escape probability routines */
108  DrvEscP( );
109  }
110 
111  else if( nMatch("HYAS",chCard) )
112  {
113  /* option to drive Jason's hydrogen transition probabilities */
114  DrvHyas();
115  }
116 
117  else if( nMatch("GAUN",chCard) )
118  {
119  /* drive gaunt factor routine */
120  dgaunt();
121  }
122 
123  else if( nMatch("POIN",chCard) )
124  {
125  /* later on, check cell pointers, centers, widths */
126  fprintf( ioQQQ, " Define entire model first, later I will ask for energies.\n" );
127  trace.lgPtrace = true;
128  }
129 
130  else if( nMatch("PUMP",chCard) )
131  {
132  char chInput[INPUT_LINE_LENGTH];
133  lgEOL = false;
134  fprintf( ioQQQ, " Continuum pump ParseDrive entered - Enter log tau\n" );
135  while( !lgEOL )
136  {
137  if( read_whole_line( chInput , (int)sizeof(chInput) , ioStdin ) == NULL )
138  {
139  fprintf( ioQQQ, " ParseDrive.dat error getting magic number\n" );
140  cdEXIT(EXIT_FAILURE);
141  }
142  /* get tau */
143  i = 1;
144  fac = FFmtRead(chInput,&i,INPUT_LINE_LENGTH,&lgEOL);
145  if( lgEOL )
146  break;
147  fac = pow(10.,fac);
148  fprintf( ioQQQ, " Tau =%11.4e\n", fac );
149  TauDummy.Emis->TauIn = (realnum)fac;
150  fac = DrvContPump(&TauDummy);
151  fprintf( ioQQQ, " ContPump =%11.4e\n", fac );
152  fprintf( ioQQQ, " Enter null to stop, or another value.\n" );
153  }
154  fprintf( ioQQQ, " ContPump ParseDrive exits.\n" );
155  }
156 
157  else if( nMatch("STAR",chCard) )
158  {
159  char chInput[INPUT_LINE_LENGTH];
160  /* get starburst abundances */
161  for( i=0; i < 40; i++ )
162  {
163  zed = ((double)i+1.)/4. + 0.01;
164  sprintf( chInput, "starburst, zed=%10.4f", zed );
165  abund_starburst(chInput);
166  fprintf( ioQQQ, "%8.1e", zed );
167  for(n=0; n < LIMELM; n++)
168  fprintf( ioQQQ, "%8.1e", abund.solar[n] );
169  fprintf( ioQQQ, "\n" );
170  }
171  }
172 
173  else
174  {
175  fprintf( ioQQQ,
176  " Unrecognized key; keys are FFmtRead, CASE, GAUNT, HYAS, POINT, MOLE, STAR, "
177  "PUMP, and ESCApe. Sorry.\n" );
178  cdEXIT(EXIT_FAILURE);
179  }
180  return;
181 }
182 
183 /*DrvEscP user queries escape probability routines, which return values */
184 STATIC void DrvEscP( void )
185 {
186  char chCard[INPUT_LINE_LENGTH];
187  bool lgEOL;
188  long i;
189  double tau;
190 
191  DEBUG_ENTRY( "DrvEscP()" );
192 
193  /* this routine is enterd with the command escape probability, and
194  * drives the escape probability routine to compare answers */
195  fprintf( ioQQQ, " Enter the log of the one-sided optical depth; line with no number to stop.\n" );
196 
197  lgEOL = false;
198  while( !lgEOL )
199  {
200  if( read_whole_line( chCard , (int)sizeof(chCard) , ioStdin ) == NULL )
201  {
202  break;
203  }
204 
205  i = 1;
206  tau = FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
207  if( lgEOL )
208  {
209  break;
210  }
211 
212  tau = pow(10.,tau);
213  fprintf( ioQQQ, "tau was %e\n", tau );
214  fprintf( ioQQQ, " ESCINC=%11.3e\n", esc_PRD_1side(tau,1e-4) );
215  fprintf( ioQQQ, " ESCCOM=%11.3e\n", esc_CRDwing_1side(tau,1e-4 ) );
216  fprintf( ioQQQ, " ESCA0K2=%11.3e\n", esca0k2(tau) );
217 
218  }
219  return;
220 }
221 
222 /*DrvCaseBHS allow user to query hydrogen A's, asks for up, low level, gives A, drive hyas */
223 STATIC void DrvCaseBHS(void)
224 {
225  char chCard[INPUT_LINE_LENGTH];
226  bool lgEOL,
227  lgHit;
228  long int i,
229  n1,
230  nelem ,
231  n2;
232  double Temperature,
233  Density;
234 
235  DEBUG_ENTRY( "DrvCaseBHS()" );
236 
237  /* this routine is entered with the command DRIVE CASEB, and
238  * interpolates on the Hummer & Storey case b data set */
239 
240  /* read in some external data files, but only if this is first call */
241  fprintf(ioQQQ," I will get needed H data files. This will take a second.\n");
242  atmdat_readin();
243 
244  {
245  /* following should be set true to print input lines */
246  /*@-redef@*/
247  enum {DEBUG_LOC=false};
248  /*@+redef@*/
249  if( DEBUG_LOC )
250  {
251  double xLyman , alpha;
252  long int ipHi;
253  nelem = 2;
254  Temperature = 2e4;
255  Density = 1e2;
256  for( ipHi=3; ipHi<25; ++ipHi )
257  {
258  double photons = (1./POW2(ipHi-1.)-1./POW2((double)ipHi) ) /(1.-1./ipHi/ipHi );
259  xLyman = atmdat_HS_caseB(1,ipHi, nelem,Temperature , Density , 'A' );
260  alpha = atmdat_HS_caseB(ipHi-1,ipHi, nelem,Temperature , Density , 'A' );
261  fprintf(ioQQQ,"%li\t%.3e\t%.3e\n", ipHi, xLyman/alpha*photons, photons );
262  }
263  cdEXIT(EXIT_SUCCESS);
264  }
265  }
266 
267  /* first get the charge, only H and He at present */
268  lgHit = false;
269  nelem = 0;
270  while( !lgHit )
271  {
272  fprintf( ioQQQ, " Enter atomic number of species, either 1(H) or 2(He).\n" );
273  if( read_whole_line( chCard , (int)sizeof(chCard) , ioStdin ) == NULL )
274  {
275  fprintf( ioQQQ, " error getting species \n" );
276  cdEXIT(EXIT_FAILURE);
277  }
278 
279  i = 1;
280  nelem = (long int)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
281  if( lgEOL || nelem< 1 || nelem > 2 )
282  {
283  fprintf( ioQQQ, " must be either 1 or 2!\n" );
284  }
285  else
286  {
287  lgHit = true;
288  }
289  }
290 
291  fprintf(ioQQQ," In the following temperatures <10 are log, >=10 linear.\n");
292  fprintf(ioQQQ," The density is always a log.\n");
293  fprintf(ioQQQ," The order of the quantum numbers do not matter.\n");
294  fprintf(ioQQQ," The smallest must not be smaller than 2,\n");
295  fprintf(ioQQQ," and the largest must not be larger than 25.\n");
296  fprintf(ioQQQ," Units of emissivity are erg cm^3 s^-1\n\n");
297  fprintf(ioQQQ," The limits of the HS tables are 2 <= n <= 25.\n");
298 
299  lgHit = true;
300  /* this is always true */
301  while( lgHit )
302  {
303  fprintf( ioQQQ, " Enter 4 numbers, temperature, density, 2 quantum numbers, null line stop.\n" );
304  if( read_whole_line( chCard , (int)sizeof(chCard) , ioStdin ) == NULL )
305  {
306  fprintf( ioQQQ, " Thanks for interpolating on the Hummer & Storey data set!\n" );
307  cdEXIT(EXIT_FAILURE);
308  }
309 
310  i = 1;
311  Temperature = FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
312  if( lgEOL )
313  {
314  fprintf( ioQQQ, " error getting temperature!\n" );
315  break;
316  }
317 
318  /* log if less than 10 */
319  if( Temperature < 10. )
320  {
321  Temperature = pow(10., Temperature );
322  }
323  fprintf(ioQQQ," Temperature is %g\n", Temperature );
324 
325  Density = FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
326  if( lgEOL )
327  {
328  fprintf( ioQQQ, " error getting density!\n" );
329  break;
330  }
331  Density = pow(10., Density );
332  fprintf(ioQQQ," Density is %g\n", Density );
333 
334  /* these quantum numbers can be in any order */
335  n1 = (long)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
336  if( lgEOL )
337  {
338  fprintf( ioQQQ, " error getting quantum number!\n" );
339  break;
340  }
341 
342  n2 = (long)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
343  if( lgEOL )
344  {
345  fprintf( ioQQQ, " error getting quantum number!\n" );
346  break;
347  }
348 
349  if( MAX2( n1 , n2 ) > 25 )
350  {
351  fprintf( ioQQQ," The limits of the HS tables are 2 <= n <= 25. Sorry.\n");
352  break;
353  }
354 
355  fprintf( ioQQQ,
356  " 4pJ(%ld,%ld)/n_e n_p=%11.3e\n",
357  n1, n2,
358  atmdat_HS_caseB(n1,n2, nelem,Temperature , Density , 'B' ) );
359 
360  /* this is check that we were in bounds */
361 #if 0
362  {
363  long j;
364  double tempTable[33] = {
365  11870.,12490.,12820.,
366  11060.,17740.,12560.,
367  16390.,16700.,11360.,
368  10240.,20740.,12030.,
369  14450.,19510.,12550.,
370  16470.,16560.,12220.,
371  15820.,12960.,10190.,
372  12960.,14060.,12560.,
373  11030.,10770.,13360.,
374  10780.,11410.,11690.,
375  12500.,13190.,21120. };
376  double edenTable[33] = {
377  10.,270.,80.,10.,70.,
378  110.,200.,10.,40.,90.,
379  340.,80.,60.,340.,30.,
380  120.,10.,50.,450.,30.,
381  180.,20.,170.,60.,20.,
382  40.,30.,20.,100.,130.,
383  10.,10.,110. };
384 
385 
386  for( j=0; j<33; j++ )
387  {
388  double halpha, hbeta, hgamma;
389 
390  halpha = atmdat_HS_caseB(2,3, 1,tempTable[j] , edenTable[j] , 'B' );
391  hbeta = atmdat_HS_caseB(2,4, 1,tempTable[j] , edenTable[j] , 'B' );
392  hgamma = atmdat_HS_caseB(2,5, 1,tempTable[j] , edenTable[j] , 'B' );
393 
394  fprintf( ioQQQ, "%e\t%e\t%e\t%e\n",
395  tempTable[j],
396  edenTable[j],
397  halpha/hbeta,
398  hgamma/hbeta );
399  }
400  }
401 #endif
402  }
403 
404  fprintf( ioQQQ, " Thanks for interpolating on the Hummer & Storey data set!\n" );
405  cdEXIT(EXIT_FAILURE);
406 
407 }
408 
409 /*DrvHyas allow user to query hydrogen A's, asks for up, low level, gives A, drive hyas */
410 STATIC void DrvHyas(void)
411 {
412  char chCard[INPUT_LINE_LENGTH];
413  bool lgEOL;
414  long int i, nHi, lHi, nLo, lLo;
415 
416  DEBUG_ENTRY( "DrvHyas()" );
417 
418  /* this routine is entered with the command DRIVE HYAS, and
419  * drives Jason's hydrogen einstein A routines */
420 
421  nHi = 1;
422  /* nHi never lt 1 */
423  while( nHi != 0 )
424  {
425  fprintf( ioQQQ, " Enter four quantum numbers (n, l, n', l'), null line to stop.\n" );
426  if( read_whole_line( chCard , (int)sizeof(chCard) , ioStdin ) == NULL )
427  {
428  fprintf( ioQQQ, " error getting drvhyas \n" );
429  cdEXIT(EXIT_FAILURE);
430  }
431 
432  i = 1;
433  nHi = (long int)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
434  if( lgEOL )
435  break;
436 
437  lHi = (long int)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
438  if( lgEOL )
439  {
440  fprintf( ioQQQ, " must be four numbers!\n" );
441  break;
442  }
443 
444  if( lHi >= nHi )
445  {
446  fprintf( ioQQQ, " l must be less than n!\n" );
447  break;
448  }
449 
450  nLo = (long int)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
451  if( lgEOL )
452  {
453  fprintf( ioQQQ, " must be four numbers!\n" );
454  break;
455  }
456 
457  lLo = (long int)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
458  if( lgEOL )
459  {
460  fprintf( ioQQQ, " must be four numbers!\n" );
461  break;
462  }
463 
464  if( lLo >= nLo )
465  {
466  fprintf( ioQQQ, " l must be less than n!\n" );
467  break;
468  }
469 
470  if( nLo > nHi )
471  {
472  long nTemp, lTemp;
473 
474  /* swap hi and lo */
475  nTemp = nLo;
476  lTemp = lLo;
477  nLo = nHi;
478  lLo = lHi;
479  nHi = nTemp;
480  lHi = lTemp;
481  }
482 
483  fprintf( ioQQQ, " A(%3ld,%3ld->%3ld,%3ld)=%11.3e\n",
484  nHi, lHi, nLo, lLo,
485  H_Einstein_A( nHi, lHi, nLo, lLo, 1 ) );
486 
487  }
488  fprintf( ioQQQ, " Driver exits, enter next line.\n" );
489 
490  return;
491 }
492 
493 /*dgaunt drive gaunt factor routines by letting user query values */
494 STATIC void dgaunt(void)
495 {
496  char chCard[INPUT_LINE_LENGTH];
497  bool lgEOL;
498  int inputflag;
499  long int i,
500  ierror;
501  realnum enerlin[1];
502  double SaveTemp;
503  double z,mygaunt=0.;
504  double loggamma2, logu;
505 
506  DEBUG_ENTRY( "dgaunt()" );
507 
508  SaveTemp = phycon.te;
509 
510  /* this routine is entered with the command DRIVE GAUNT, and
511  * drives the gaunt factor routine to check range
512  * */
513  fprintf( ioQQQ, " Enter 0 to input temp, energy, and net charge, or 1 for u and gamma**2.\n" );
514  if( read_whole_line( chCard , (int)sizeof(chCard) , ioStdin ) == NULL )
515  {
516  fprintf( ioQQQ, " dgaunt error getting magic number\n" );
517  cdEXIT(EXIT_FAILURE);
518  }
519  i = 1;
520  inputflag = (int)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
521 
522  if( inputflag == 0 )
523  {
524  fprintf( ioQQQ, " Enter the temperature (log if <=10), energy (Ryd), and net charge. Null line to stop.\n" );
525  /* >>chng 96 july 07, got rid of statement labels replacing with do while
526  * */
527  ierror = 0;
528  while( ierror == 0 )
529  {
530  if( read_whole_line( chCard , (int)sizeof(chCard) , ioStdin ) == NULL )
531  {
532  fprintf( ioQQQ, " dgaunt error getting magic number\n" );
533  cdEXIT(EXIT_FAILURE);
534  }
535  i = 1;
536  phycon.alogte = FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
537  /* the line may be trash but ierror will pick it up */
538  if( lgEOL )
539  {
540  fprintf( ioQQQ, " Gaunt driver exits, enter next line.\n" );
541  break;
542  }
543  /* numbers less than or equal to 10 are the log of the temperature */
544  double TeNew;
545  if( phycon.alogte > 10. )
546  {
547  TeNew = phycon.alogte;
548  }
549  else
550  {
551  TeNew = pow(10.,phycon.alogte);
552  }
553  TempChange(TeNew , false);
554 
555  enerlin[0] = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
556  if( lgEOL || enerlin[0] == 0. )
557  {
558  fprintf( ioQQQ, " Sorry, but there should be two more numbers, energy and charge.\n" );
559  }
560 
561  z = (double)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
562  if( lgEOL || z == 0. )
563  {
564  fprintf( ioQQQ, " Sorry, but there should be a third number, charge.\n" );
565  }
566 
567  /* This is non-thermally averaged gaunt factors. */
568  mygaunt = cont_gaunt_calc( (double)phycon.te, z, enerlin[0] );
569 
570  fprintf( ioQQQ, " Using my routine, Gff= \t" );
571  fprintf( ioQQQ, "%11.3e\n", mygaunt );
572 
573  }
574  }
575  else
576  {
577  /* this routine is entered with the command DRIVE GAUNT, and
578  * drives the gaunt factor routine to check range
579  * */
580  fprintf( ioQQQ, " Enter log u and log gamma2. Null line to stop.\n" );
581  ierror = 0;
582  while( ierror == 0 )
583  {
584  if( read_whole_line( chCard , (int)sizeof(chCard) , ioStdin ) == NULL )
585  {
586  fprintf( ioQQQ, " dgaunt error getting magic number\n" );
587  cdEXIT(EXIT_FAILURE);
588  }
589  i = 1;
590  logu = FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
591  /* the line may be trash but ierror will pick it up */
592  if( lgEOL )
593  {
594  fprintf( ioQQQ, " Gaunt driver exits, enter next line.\n" );
595  break;
596  }
597 
598  loggamma2 = FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
599  if( lgEOL )
600  {
601  fprintf( ioQQQ, " Sorry, but there should be another numbers, log gamma2.\n" );
602  }
603 
604  /* This is my attempt to calculate non-thermally averaged gaunt factors. */
605  mygaunt = cont_gaunt_calc( TE1RYD/pow(10.,loggamma2), 1., pow(10.,logu-loggamma2) );
606 
607  TempChange(TE1RYD/pow(10.,loggamma2) , false);
608 
609  fprintf( ioQQQ, " Using my routine, Gaunt factor is" );
610  fprintf( ioQQQ, "%11.3e\n", mygaunt );
611  }
612  }
613 
614  TempChange(SaveTemp , false);
615  return;
616 }

Generated for cloudy by doxygen 1.8.1.1