cloudy  trunk
 All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
mole_co_etc.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 /*CO_Init called from cdInit to initialize co routines */
4 /*CO_update_chem_rates update rate coefficients, only temp part - in mole_co_etc.c */
5 #include "cddefines.h"
6 #include "physconst.h"
7 #include "mole.h"
8 #include "mole_co_priv.h"
9 #include "hmi.h"
10 #include "rfield.h"
11 #include "dense.h"
12 #include "ionbal.h"
13 #include "grainvar.h"
14 #include "timesc.h"
15 /*lint -e778 constant expression evaluates to 0 in operation '-' */
16 
17 /*CO_update_chem_rates update rate coefficients, only temp part - in mole_co_etc.c
18  * called in conv_base before any chemistry or ionization is done */
19 
22 
23 STATIC void newelement(const char label[], int ipion,
24  int priority);
25 STATIC struct molecule *newspecies(const char label[7], enum spectype type,
26  enum molstate state, realnum *location, double frac0);
27 STATIC struct chem_element_s *findelement(const char buf[]);
28 STATIC int isactive(data_u *dat);
29 STATIC int ispassive(data_u *dat);
30 STATIC int isCOnet(data_u *dat);
31 
33 /* List of element structures indexed by atom index
34  -- could use (element_list[nelem] != NULL) for mole.lgElem_in_chemistry[nelem] */
36 int32 *ipiv;
38 
41 
42 /*=================================================================*/
43 /*CO_Init called from cdInit to initialize CO routines */
44 void CO_Init(void)
45 {
46  long int i,
47  nelem;
48  struct molecule *sp;
49  static bool lgCO_Init_called=false;
50 
51  DEBUG_ENTRY( "CO_Init()" );
52 
53  /* these tell the molecular solver what zone and iteration it has
54  * been evaluated on */
55  co.co_nzone = -2;
56  co.iteration_co = -2;
57 
58  /* prevent memory leaks */
59  /* \todo this is a temporary fix for PR14. We should improve the overall design
60  * of this code to prevent valid pointers being overwritten in a second call to CO_Init */
61  if( lgCO_Init_called )
62  {
63  return;
64  }
65 
66  /* say that we have been called */
67  lgCO_Init_called = true;
68 
69  mole_priv.spectab = newhash(NULL);
70  mole_priv.reactab = newhash(NULL);
71  mole_priv.elemtab = newhash(NULL);
72 
73  for(nelem=ipHYDROGEN; nelem<LIMELM; ++nelem )
74  {
75  element_list[nelem] = NULL;
76  mole.lgElem_in_chemistry[nelem] = false;
77  }
78 
79  /* set up concordance of elemental species to external Cloudy indices */
80  /* final number is 'element priority':-
81  rjrw 2006 Aug 11: Nick Abel explains this order is roughly gas phase
82  abundance -- the encoding of atomic species is from the last
83  elements of "enum molecule_codes" in mole.h. Should this make any
84  difference? In practice it does. */
85  newelement("C" ,ipCARBON,3);
86  newelement("O" ,ipOXYGEN,2);
87  newelement("Si",ipSILICON,5);
88  newelement("N" ,ipNITROGEN,4);
89  newelement("S" ,ipSULPHUR,6);
90  newelement("Cl",ipCHLORINE,7);
91  newelement("H" ,ipHYDROGEN,1);
92  newelement("He",ipHELIUM,0);
93  newelement("Mg",ipMAGNESIUM,0);
94  newelement("Fe",ipIRON,0);
95 
96  /* set up properties of molecular species -- chemical formulae,
97  array indices, elementary components (parsed from formula),
98  status within CO network, location of stored value external
99  to CO network (must be floating point). */
100 
101  /* final arguments are relative to neutral carbon and are the solution to the first zone of the
102  * TH85 pdr */
103 
104  /* Sizes of different parts of network are calculated in newspecies */
106  null_mole.hevmol = null_mole.hevcol = 0.; /* Non-molecule to allow valid pointer return inactive species are accessed */
107  null_mole.index = -1;
108 
109  sp = newspecies("C ",MOLECULE,ACTIVE,NULL,1.00e+00);
110  sp = newspecies("O ",MOLECULE,ACTIVE,NULL,1.13e+05);
111  sp = newspecies("Si ",MOLECULE,ACTIVE,NULL,3.56e-04);
112  sp = newspecies("N ",MOLECULE,ACTIVE,NULL,2.25e+00);
113  sp = newspecies("S ",MOLECULE,ACTIVE,NULL,6.90e-03);
114  sp = newspecies("Cl ",MOLECULE,ACTIVE,NULL,6.90e-03);
115  sp = newspecies("C+ ",MOLECULE,ACTIVE,NULL,6.79e+04);
116  sp = newspecies("O+ ",MOLECULE,ACTIVE,NULL,1.58e+01);
117  sp = newspecies("Si+ ",MOLECULE,ACTIVE,NULL,1.79e+02);
118  sp = newspecies("N+ ",MOLECULE,ACTIVE,NULL,2.62e-10);
119  sp = newspecies("S+ ",MOLECULE,ACTIVE,NULL,1.79e+03);
120  sp = newspecies("Cl+ ",MOLECULE,ACTIVE,NULL,3.71e-09);
121  sp = newspecies("CH ",MOLECULE,ACTIVE,NULL,1.10e-06);
122  sp = newspecies("CH+ ",MOLECULE,ACTIVE,NULL,6.99e-05);
123  sp = newspecies("OH ",MOLECULE,ACTIVE,NULL,2.15e-03);
124  sp = newspecies("OH+ ",MOLECULE,ACTIVE,NULL,2.45e-04);
125  sp = newspecies("O2 ",MOLECULE,ACTIVE,NULL,1.91e-07);
126  sp = newspecies("CO ",MOLECULE,ACTIVE,NULL,4.05e-05);
127  sp = newspecies("CO+ ",MOLECULE,ACTIVE,NULL,2.16e-06);
128  sp = newspecies("H2O ",MOLECULE,ACTIVE,NULL,1.60e-11);
129  sp = newspecies("H2O+ ",MOLECULE,ACTIVE,NULL,1.27e-10);
130  sp = newspecies("O2+ ",MOLECULE,ACTIVE,NULL,1.17e-06);
131  sp = newspecies("H3O+ ",MOLECULE,ACTIVE,NULL,3.44e-17);
132  sp = newspecies("CH2+ ",MOLECULE,ACTIVE,NULL,3.71e-09);
133  sp = newspecies("CH2 ",MOLECULE,ACTIVE,NULL,8.59e-13);
134  sp = newspecies("HCO+ ",MOLECULE,ACTIVE,NULL,4.01e-10);
135  sp = newspecies("CH3+ ",MOLECULE,ACTIVE,NULL,7.02e-16);
136  sp = newspecies("CH3 ",MOLECULE,ACTIVE,NULL,4.59e-19);
137  sp = newspecies("CH4 ",MOLECULE,ACTIVE,NULL,9.12e-28);
138  sp = newspecies("CH4+ ",MOLECULE,ACTIVE,NULL,1.03e-28);
139  sp = newspecies("CH5+ ",MOLECULE,ACTIVE,NULL,8.16e-28);
140  sp = newspecies("SiH2+ ",MOLECULE,ACTIVE,NULL,2.07e-13);
141  sp = newspecies("SiH ",MOLECULE,ACTIVE,NULL,1.80e-14);
142  sp = newspecies("SiOH+ ",MOLECULE,ACTIVE,NULL,5.39e-15);
143  sp = newspecies("SiO ",MOLECULE,ACTIVE,NULL,3.89e-14);
144  sp = newspecies("SiO+ ",MOLECULE,ACTIVE,NULL,2.27e-08);
145  sp = newspecies("N2 ",MOLECULE,ACTIVE,NULL,2.46e-17);
146  sp = newspecies("N2+ ",MOLECULE,ACTIVE,NULL,2.22e-17);
147  sp = newspecies("NO ",MOLECULE,ACTIVE,NULL,5.99e-12);
148  sp = newspecies("NO+ ",MOLECULE,ACTIVE,NULL,1.31e-11);
149  sp = newspecies("S2 ",MOLECULE,ACTIVE,NULL,1.21e-11);
150  sp = newspecies("S2+ ",MOLECULE,ACTIVE,NULL,1.00e-13);
151  sp = newspecies("OCN ",MOLECULE,ACTIVE,NULL,3.86e-11);
152  sp = newspecies("OCN+ ",MOLECULE,ACTIVE,NULL,1.65e-22);
153  sp = newspecies("NH ",MOLECULE,ACTIVE,NULL,1.30e-09);
154  sp = newspecies("NH+ ",MOLECULE,ACTIVE,NULL,2.20e-10);
155  sp = newspecies("NH2 ",MOLECULE,ACTIVE,NULL,6.78e-20);
156  sp = newspecies("NH2+ ",MOLECULE,ACTIVE,NULL,1.71e-16);
157  sp = newspecies("NH3 ",MOLECULE,ACTIVE,NULL,1.25e-29);
158  sp = newspecies("NH3+ ",MOLECULE,ACTIVE,NULL,3.03e-23);
159  sp = newspecies("NH4+ ",MOLECULE,ACTIVE,NULL,1.15e-33);
160  sp = newspecies("CN ",MOLECULE,ACTIVE,NULL,3.38e-11);
161  sp = newspecies("CN+ ",MOLECULE,ACTIVE,NULL,5.07e-11);
162  sp = newspecies("HCN ",MOLECULE,ACTIVE,NULL,1.07e-17);
163  sp = newspecies("HCN+ ",MOLECULE,ACTIVE,NULL,9.89e-17);
164  sp = newspecies("HNO ",MOLECULE,ACTIVE,NULL,9.09e-21);
165  sp = newspecies("HNO+ ",MOLECULE,ACTIVE,NULL,1.91e-18);
166  sp = newspecies("HS ",MOLECULE,ACTIVE,NULL,1.71e-14);
167  sp = newspecies("HS+ ",MOLECULE,ACTIVE,NULL,4.83e-13);
168  sp = newspecies("CS ",MOLECULE,ACTIVE,NULL,6.69e-18);
169  sp = newspecies("CS+ ",MOLECULE,ACTIVE,NULL,4.12e-11);
170  sp = newspecies("NO2 ",MOLECULE,ACTIVE,NULL,2.67e-26);
171  sp = newspecies("NO2+ ",MOLECULE,ACTIVE,NULL,2.41e-23);
172  sp = newspecies("NS ",MOLECULE,ACTIVE,NULL,3.12e-11);
173  sp = newspecies("NS+ ",MOLECULE,ACTIVE,NULL,6.81e-13);
174  sp = newspecies("SO ",MOLECULE,ACTIVE,NULL,4.00e-15);
175  sp = newspecies("SO+ ",MOLECULE,ACTIVE,NULL,1.20e-07);
176  sp = newspecies("SiN ",MOLECULE,ACTIVE,NULL,7.03e-12);
177  sp = newspecies("SiN+ ",MOLECULE,ACTIVE,NULL,7.60e-14);
178  sp = newspecies("N2O ",MOLECULE,ACTIVE,NULL,1.05e-11);
179  sp = newspecies("HCS+ ",MOLECULE,ACTIVE,NULL,8.91e-17);
180  sp = newspecies("OCS ",MOLECULE,ACTIVE,NULL,8.83e-24);
181  sp = newspecies("OCS+ ",MOLECULE,ACTIVE,NULL,1.72e-21);
182  sp = newspecies("HNC ",MOLECULE,ACTIVE,NULL,4.00e-15);
183  sp = newspecies("HCNH+ ",MOLECULE,ACTIVE,NULL,4.00e-15);
184  sp = newspecies("C2 ",MOLECULE,ACTIVE,NULL,3.71e-09);
185  sp = newspecies("C2+ ",MOLECULE,ACTIVE,NULL,8.59e-13);
186  sp = newspecies("CCl ",MOLECULE,ACTIVE,NULL,4.00e-15);
187  sp = newspecies("ClO ",MOLECULE,ACTIVE,NULL,4.00e-15);
188  sp = newspecies("HCl+ ",MOLECULE,ACTIVE,NULL,4.00e-15);
189  sp = newspecies("HCl ",MOLECULE,ACTIVE,NULL,4.00e-15);
190  sp = newspecies("H2Cl+ ",MOLECULE,ACTIVE,NULL,4.00e-15);
191  sp = newspecies("CCl+ ",MOLECULE,ACTIVE,NULL,4.00e-15);
192  sp = newspecies("H2CCl+",MOLECULE,ACTIVE,NULL,0.00e+00);
193  sp = newspecies("ClO+ ",MOLECULE,ACTIVE,NULL,4.00e-15);
194  /* fprintf(stderr,"ETC: %d %d\n",gv.lgDustOn, mole.lgGrain_mole_deplete); */
196  {
197  sp = newspecies("COgrn ",MOLECULE,ACTIVE,NULL,1.00e-15);
198  sp = newspecies("H2Ogrn",MOLECULE,ACTIVE,NULL,0.00e+00);
199  sp = newspecies("OHgrn ",MOLECULE,ACTIVE,NULL,1.00e-15);
200  }
201  sp = newspecies("C2H ",MOLECULE,ACTIVE,NULL,4.00e-15);
202  sp = newspecies("C2H+ ",MOLECULE,ACTIVE,NULL,4.00e-15);
203  sp = newspecies("C3 ",MOLECULE,ACTIVE,NULL,1.00e-15);
204  sp = newspecies("C3+ ",MOLECULE,ACTIVE,NULL,4.00e-15);
205  sp = newspecies("C3H+ ",MOLECULE,ACTIVE,NULL,4.00e-15);
206  sp = newspecies("C3H ",MOLECULE,ACTIVE,NULL,4.00e-15);
207  sp = newspecies("C2H2 ",MOLECULE,ACTIVE,NULL,4.00e-15);
208  sp = newspecies("C2H2+ ",MOLECULE,ACTIVE,NULL,4.00e-15);
209  sp = newspecies("C2H3+ ",MOLECULE,ACTIVE,NULL,4.00e-15);
210  sp = newspecies("N2H+ ",MOLECULE,ACTIVE,NULL,4.00e-15);
211  /* Add passive species to complete network */
212  sp = newspecies("e- ",OTHER,PASSIVE,&(dense.eden_f),0.00e+00);
213  sp->nElec = -1; sp->mole_mass = (realnum)ELECTRON_MASS; /* Augment properties for this non-molecular species */
214 
215  sp = newspecies("Fe ",MOLECULE,PASSIVE,&(dense.xIonDense[ipIRON][0]),0.00e+00);
216  sp = newspecies("Fe+ ",MOLECULE,PASSIVE,&(dense.xIonDense[ipIRON][1]),0.00e+00);
217  sp = newspecies("H ",MOLECULE,PASSIVE,&(dense.xIonDense[ipHYDROGEN][0]),0.00e+00);
218  sp = newspecies("H- ",MOLECULE,PASSIVE,&(hmi.Hmolec[ipMHm]),0.00e+00);
219  sp = newspecies("H+ ",MOLECULE,PASSIVE,&(dense.xIonDense[ipHYDROGEN][1]),0.00e+00);
220  sp = newspecies("H2 ",MOLECULE,PASSIVE,&(hmi.Hmolec[ipMH2g]),0.00e+00);
221  sp = newspecies("H2* ",MOLECULE,PASSIVE,&(hmi.Hmolec[ipMH2s]),0.00e+00);
222  sp = newspecies("H2+ ",MOLECULE,PASSIVE,&(hmi.Hmolec[ipMH2p]),0.00e+00);
223  sp = newspecies("H3+ ",MOLECULE,PASSIVE,&(hmi.Hmolec[ipMH3p]),0.00e+00);
224  sp = newspecies("He ",MOLECULE,PASSIVE,&(dense.xIonDense[ipHELIUM][0]),0.00e+00);
225  sp = newspecies("He+ ",MOLECULE,PASSIVE,&(dense.xIonDense[ipHELIUM][1]),0.00e+00);
226  sp = newspecies("Mg ",MOLECULE,PASSIVE,&(dense.xIonDense[ipMAGNESIUM][0]),0.00e+00);
227  sp = newspecies("Mg+ ",MOLECULE,PASSIVE,&(dense.xIonDense[ipMAGNESIUM][1]),0.00e+00);
228 
229  /* Create linear list of species and populate it... */
230  COmole = (struct molecule **)MALLOC((size_t)mole.num_comole_tot*
231  sizeof(struct molecule *));
232 
233  /* ...first active species */
234  i = makeplist(mole_priv.spectab,(void **)COmole,
236  ASSERT (i == mole.num_comole_calc);
237 
238  /* ...then passive species at end of list */
242 
243  /* Set molecule indices to order of list just created */
244  for(i=0;i<mole.num_comole_tot;i++)
245  {
246  COmole[i]->index = i;
247  }
248 
249  /* Register the atomic ladders which are active in this network */
250  for(i=0;i<mole.num_comole_calc;i++)
251  {
252  sp = COmole[i];
253  if(sp->active && sp->n_nuclei == 1)
254  {
255  if(sp->nElec == 0)
256  {
257  element_list[sp->nelem_hevmol]->ipMl = i;
258  }
259  else if(sp->nElec == 1)
260  {
261  element_list[sp->nelem_hevmol]->ipMlP = i;
262  }
263  }
264  }
265 
266  /* Create and fill cache for looping on active atomic species */
267  chem_element = (struct chem_element_s **)
268  MALLOC((size_t)(mole.num_elements)*sizeof(struct chem_element_s *));
269  i = makeplist(mole_priv.elemtab,(void **)chem_element,mole.num_elements,isCOnet);
270  ASSERT(i == mole.num_elements);
271 
272 
273  /* Create other workspace arrays which have sizes given by the species numbers */
274  mole.amat = (double **)MALLOC((size_t)mole.num_comole_calc*sizeof(double *));
275  mole.amat[0] = (double *)MALLOC((size_t)mole.num_comole_calc*
276  mole.num_comole_calc*sizeof(double));
277  for(i=1;i<mole.num_comole_calc;i++)
278  {
280  }
281  mole.b = (double *)MALLOC((size_t)mole.num_comole_calc*sizeof(double));
282  mole.c = (double **)MALLOC((size_t)mole.num_comole_tot*sizeof(double *));
283  mole.c[0] = (double *)MALLOC((size_t)mole.num_comole_tot*
284  mole.num_comole_calc*sizeof(double));
285  for(i=1;i<mole.num_comole_tot;i++)
286  {
287  mole.c[i] = mole.c[i-1]+mole.num_comole_calc;
288  }
289  ipiv = (int32 *)MALLOC((size_t)mole.num_comole_calc*sizeof(int32));
290 
291  tot_ion = (realnum *)MALLOC((size_t)mole.num_elements*sizeof(realnum));
292 
293  return;
294 
295 }
296 /*lint +e778 constant expression evaluates to 0 in operation '-' */
297 
298 /* Fill element linking structure */
299 STATIC void newelement(const char label[], int ipion, int priority)
300 {
301  struct chem_element_s *element;
302  data_u *p;
303  int exists;
304 
305  DEBUG_ENTRY("newelement()");
306 
307  element = (struct chem_element_s *) MALLOC(sizeof(struct chem_element_s));
308  element->ipCl = ipion;
309  element->ipZ = priority;
310  ASSERT ((size_t)strlen(label) < sizeof(element->chName));
311  strncpy(element->chName,label,sizeof(element->chName));
312  element->ipMl = element->ipMlP = -1; /* Chemical network species indices not yet defined */
313  p = addentry(element->chName,0,mole_priv.elemtab,&exists);
314  p->p = (void *) element;
315  element_list[ipion] = element;
316 }
317 /* Parse species string to find constituent atoms, charge etc. */
318 #include <string.h>
319 #include <ctype.h>
320 STATIC struct molecule *newspecies(const char label[7], enum spectype type,
321  enum molstate state, realnum *location, double frac0)
322 {
323  int exists;
324  data_u *p;
325  char mylab[7], thisel[3], *s;
326  long int i, n, nnuc, nelem, nel;
327  struct molecule *mol;
328  struct chem_element_s *el, *maxel;
329 
330  DEBUG_ENTRY("newspecies()");
331 
332  mol = (struct molecule *) MALLOC(sizeof(struct molecule));
333 
335  if(state == ACTIVE)
337 
338  mol->pdr_mole_co = (realnum) frac0;
339  for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem )
340  {
341  mol->nElem[nelem] = 0;
342  }
343  mol->co_save = 0.;
344 
345  strncpy(mol->label,label,sizeof(mol->label));
346  s = strchr(mol->label,' ');
347  if(s)
348  *s = '\0';
349 
350  /* Insert species into hash table */
351  p = addentry(mol->label,0,mole_priv.spectab,&exists);
352  p->p = (void *) mol;
353 
354  if(state == ACTIVE)
355  {
356  mol->active = 1;
357  }
358  else
359  {
360  mol->active = 0;
361  }
362  mol->location = location;
363  mol->n_nuclei = 0;
364 
365  if(type == MOLECULE)
366  {
367  /* Copy to private workspace */
368  strncpy(mylab,label,7);
369 
370  /* Trailing modifiers */
371  mol->nElec = mol->Excit = 0;
372 
373  /* Excitation... */
374  s = strpbrk(mylab,"*");
375  if(s)
376  {
377  mol->Excit = 1;
378  *s = '\0';
379  }
380 
381  /* ...Charge */
382  s = strpbrk(mylab,"+-");
383  if(s)
384  {
385  if(isdigit(*(s+1)))
386  n = atoi(s+1);
387  else
388  n = 1;
389  if(*s == '+')
390  mol->nElec = n;
391  else
392  mol->nElec = -n;
393  *s = '\0';
394  }
395  /* ...Grain */
396  s = strstr(mylab,"grn");
397  if(s)
398  {
399  mol->lgGas_Phase = false;
400  *s = '\0';
401  }
402  else
403  {
404  mol->lgGas_Phase = true;
405  }
406 
407  /* Now analyse chemical formula */
408  nnuc = 0; /* Keep account of number of atoms contained */
409  i = 0;
410  maxel = NULL;
411  while (mylab[i] != '\0' && mylab[i] != ' ' && mylab[i] != '*')
412  {
413  /* Select next element in species, matches regexp [A-Z][a-z]? */
414  nel = 0;
415  thisel[nel++] = mylab[i++];
416  if(islower(mylab[i]))
417  {
418  thisel[nel++] = mylab[i++];
419  }
420  thisel[nel] = '\0';
421 
422  el = findelement(thisel);
423  if(el == NULL)
424  {
425  fprintf(stderr,"Did not recognize element at %s in \"%s \"[%ld]\n",mylab+i,label,i);
426  cdEXIT( EXIT_FAILURE );
427  }
428 
429  /* Determine 'heaviest' atom in molecular species */
430  /* >> chng 06 Sep 25 rjrw: only C, O, Si, N, S, Cl and H count for this */
431  if(el->ipZ != 0 && (maxel == NULL || el->ipZ > maxel->ipZ))
432  {
433  maxel = el;
434  }
435 
436  if(isdigit(mylab[i])) /* If there is >1 of this atom */
437  {
438  n = 0;
439  do {
440  n = 10*n+(long int)(mylab[i]-'0');
441  i++;
442  } while (isdigit(mylab[i]));
443  }
444  else
445  {
446  n = 1;
447  }
448  mol->nElem[el->ipCl] += n;
449  nnuc += n;
450  }
451 
452  mol->n_nuclei = nnuc;
453 
454  ASSERT(mol->active == 0 || maxel != NULL );
455 
456  if(maxel != NULL)
457  mol->nelem_hevmol = maxel->ipCl;
458  else
459  mol->nelem_hevmol = -1;
460 
461  mol->mole_mass = 0.;
462  for(i=0;i<LIMELM;++i)
463  {
464  mol->mole_mass +=
466  }
467 
468  /* If we've found a neutral atomic species active in this network */
469  if(mol->n_nuclei == 1 && mol->active && mol->nElec == 0)
470  {
471  mole.num_elements++;
473  }
474  }
475  return mol;
476 }
478 {
479 
480  DEBUG_ENTRY("isactive()");
481 
482  struct molecule *p = (struct molecule *) (dat->p);
483 
484  return p->active;
485 }
487 {
488  return !((struct molecule *) dat->p)->active;
489 }
491 {
492  return (((struct chem_element_s *) dat->p)->ipMl != -1);
493 }
494 
495 struct molecule *findspecies(const char buf[])
496 {
497  data_u *p;
498 
499  DEBUG_ENTRY("findspecies()");
500 
501  p = lookup(buf,0,mole_priv.spectab);
502 
503  if(p)
504  return (struct molecule *) p->p;
505  else
506  return &null_mole;
507 }
508 
509 STATIC struct chem_element_s *findelement(const char buf[])
510 {
511  data_u *p;
512 
513  DEBUG_ENTRY("findelement()");
514 
515  p = lookup(buf,0,mole_priv.elemtab);
516 
517  if(p)
518  return (struct chem_element_s *) p->p;
519  else
520  return NULL;
521 }
522 
523 struct COmole_rate_s *CO_findrate_s(const char buf[])
524 {
525  data_u *p;
526 
527  DEBUG_ENTRY("CO_findrate_s()");
528 
529  p = lookup(buf,0,mole_priv.reactab);
530 
531  if(p)
532  return (struct COmole_rate_s *) p->p;
533  else
534  return NULL;
535 }
536 double CO_findrk(const char buf[])
537 {
538  struct COmole_rate_s *rate;
539 
540  DEBUG_ENTRY("CO_findrk()");
541 
542  rate = CO_findrate_s(buf);
543 
544  if(!rate)
545  return 0.;
546 
547  /* check for NaN */
548  ASSERT( !isnan( rate->rk ) );
549 
550  return rate->rk;
551 }
552 double CO_findrate(const char buf[])
553 {
554  struct COmole_rate_s *rate;
555  double ret;
556  int i;
557 
558  DEBUG_ENTRY("CO_findrate()");
559 
560  rate = CO_findrate_s(buf);
561  if(!rate)
562  {
563  return 0.;
564  }
565 
566  ret = rate->rk;
567  for(i=0;i<rate->nrates;i++)
568  ret *= rate->rate_species[i]->hevmol;
569  return ret;
570 }
571 
573 {
574  int i;
575 
576  DEBUG_ENTRY("CO_update_species_cache()");
577 
578  dense.eden_f = (realnum)dense.eden; /* Need floating point version for compatibility with all other values */
579 
580  for(i=0;i<mole.num_comole_tot;i++)
581  {
582  if(COmole[i]->location)
583  {
584  COmole[i]->hevmol = *(COmole[i]->location);
585 
586  /* check for NaN */
587  ASSERT( !isnan( COmole[i]->hevmol ) );
588  }
589  }
590 }
591 
592 /* Calculate rate at which molecular network abstracts species */
593 
594 /* Need to check rate_species vs reactant behaviour */
595 double CO_sink_rate(const char chSpecies[7])
596 {
597  int ipthis, i, n, nt;
598  double ratev, ratevi;
599  struct COmole_rate_s *rate;
600  struct molecule *sp;
601 
602  DEBUG_ENTRY("CO_sink_rate()");
603 
604  sp = findspecies(chSpecies);
605  ratev = 0;
606  nt = coreactions.n;
607 
608  for(n=0;n<nt;n++)
609  {
610  rate = coreactions.list[n];
611  ipthis = -1;
612  for(i=0;i<rate->nrates && ipthis == -1;i++)
613  {
614  if(rate->rate_species[i] == sp)
615  {
616  ipthis = i;
617  }
618  }
619  if(ipthis != -1) {
620  ratevi = rate->rk;
621  for(i=0;i<rate->nrates;i++)
622  {
623  if(i!=ipthis)
624  {
625  ratevi *= rate->rate_species[i]->hevmol;
626  }
627  }
628  ratev += ratevi;
629  }
630  }
631 
632  return ratev;
633 }
634 #ifdef PRINTREACT
635 STATIC void printreact(struct COmole_rate_s *rate)
636 {
637  int i;
638 
639  DEBUG_ENTRY("printreact()");
640 
641  for(i=0;i<rate->nreactants;i++) {
642  fprintf(stderr,"%s,",rate->reactants[i]->label);
643  }
644  fprintf(stderr,"=>");
645  for(i=0;i<rate->nproducts;i++) {
646  fprintf(stderr,"%s,",rate->products[i]->label);
647  }
648  fprintf(stderr,"\n");
649 
650 }
651 #endif
652 void CO_update_rks(void)
653 {
654  int n, nt;
655  struct COmole_rate_s *rate;
656 
657  DEBUG_ENTRY("CO_update_rks()");
658 
659  nt = coreactions.n;
660 
661  for(n=0;n<nt;n++)
662  {
663  rate = coreactions.list[n];
664  if(rate->fun != NULL) {
665  rate->rk = rate->a*rate->fun(rate);
666  }
667  }
668 }
669 
670 double CO_dissoc_rate(const char chSpecies[7])
671 {
672  int ipthis, i, n, nt;
673  double ratev, ratevi;
674  struct COmole_rate_s *rate;
675  struct molecule *sp;
676 
677  DEBUG_ENTRY("CO_dissoc_rate()");
678 
679  sp = findspecies(chSpecies);
680  ratev = 0;
681  nt = coreactions.n;
682 
683  for(n=0;n<nt;n++)
684  {
685  rate = coreactions.list[n];
686  if(rate->photon == -1)
687  {
688  ipthis = 0;
689  for(i=0;i<rate->nproducts;i++)
690  {
691  if(rate->products[i] == sp)
692  {
693  ipthis++;
694  }
695  }
696  if(ipthis)
697  {
698  ratevi = rate->rk;
699  for(i=0;i<rate->nrates;i++)
700  {
701  ratevi *= rate->rate_species[i]->hevmol;
702  }
703  ratev += ipthis*ratevi;
704  }
705  }
706  }
707 
708  return ratev;
709 }
710 double CO_source_rate(const char chSpecies[7])
711 {
712  int ipthis, i, n, nt;
713  double ratev, ratevi;
714  struct COmole_rate_s *rate;
715  struct molecule *sp;
716 
717  DEBUG_ENTRY("CO_source_rate()");
718 
719  sp = findspecies(chSpecies);
720  ratev = 0;
721  nt = coreactions.n;
722 
723  for(n=0;n<nt;n++)
724  {
725  rate = coreactions.list[n];
726  ipthis = 0;
727  for(i=0;i<rate->nproducts;i++)
728  {
729  if(rate->products[i] == sp)
730  {
731  ipthis++;
732  }
733  }
734  if(ipthis) {
735  ratevi = rate->rk;
736  for(i=0;i<rate->nrates;i++)
737  {
738  ratevi *= rate->rate_species[i]->hevmol;
739  }
740  ratev += ipthis*ratevi;
741  }
742  }
743 
744  return ratev;
745 }
746 
747 /* Punch all rates involving specified species */
748 void CO_punch_mol(FILE *punit, const char chSpecies[], char header[], double depth)
749 {
750  int n, nt, i, ipthis;
751  double ratevi;
752  struct COmole_rate_s *rate;
753  struct molecule *sp;
754  char *s;
755 
756  DEBUG_ENTRY("CO_punch_mol()");
757 
758  s = header;
759 
760  sp = findspecies(chSpecies);
761  if(punit == NULL)
762  {
763  sprintf (s,"#Depth");
764  s += strlen(s);
765  }
766  else
767  {
768  fprintf (punit,"%.5e",depth);
769  }
770 
771  nt = coreactions.n;
772  for(n=0;n<nt;n++) {
773  rate = coreactions.list[n];
774  ipthis = 0;
775  for(i=0;i<rate->nrates;i++)
776  {
777  if(rate->rate_species[i] == sp)
778  {
779  ipthis++;
780  }
781  }
782  for(i=0;i<rate->nproducts;i++)
783  {
784  if(rate->products[i] == sp)
785  {
786  ipthis++;
787  }
788  }
789  if(ipthis)
790  {
791  if(punit == NULL)
792  {
793  sprintf(s,"\t%s",rate->label);
794  s += strlen(s);
795  }
796  else
797  {
798  ratevi = rate->rk;
799  for(i=0;i<rate->nrates;i++)
800  {
801  ratevi *= rate->rate_species[i]->hevmol;
802  }
803  fprintf(punit,"\t%.3e",ratevi);
804  }
805  }
806  }
807  if(punit == NULL)
808  {
809  sprintf(s,"\n");
810  }
811  else
812  {
813  fprintf(punit,"\n");
814  }
815 }
816 
817 void CO_zero(void)
818 {
819  long int i;
820 
821  static bool lgFirstCall = true;
822  static long int num_comole_calc_MALLOC=-1;
823 
824  DEBUG_ENTRY("CO_zero()");
825 
826  if( lgFirstCall )
827  {
828  /* do one-time malloc of timesc.AgeCOMoleDest */
829  timesc.AgeCOMoleDest = (double*)MALLOC( mole.num_comole_calc*sizeof(double) );
830 
831  lgFirstCall = false;
832  num_comole_calc_MALLOC = mole.num_comole_calc;
833  }
834  else if( mole.num_comole_calc>num_comole_calc_MALLOC )
835  {
836  /* number of species has increased since last time - this can't happen
837  * tsuite / programs / comp4 has 95 first time, 98 second time */
838  fprintf(ioQQQ,"DISASTER - the number of species in the CO network has increased. This is not allowed.\n");
839  fprintf(ioQQQ,"This could happen if an element was initially turned off or grains not included, then the element or grains was included. There are not allowed.\n");
840  fprintf(ioQQQ,"Sorry.\n");
841  cdEXIT(EXIT_FAILURE);
842  }
843 
844  for( i=0; i < mole.num_comole_calc; i++ )
845  {
846  timesc.AgeCOMoleDest[i] = 0.;
847  }
848 
849  /* largest fraction of atoms in molecules */
850  for( i=0; i<mole.num_comole_calc; ++i )
851  COmole[i]->xMoleFracMax = 0.;
852 
853  /* zero out molecular species */
854  for( i=0; i < mole.num_comole_tot; i++ )
855  {
856  COmole[i]->hevmol = 0.;
857  COmole[i]->hevcol = 0.;
858  }
859 }

Generated for cloudy by doxygen 1.8.3.1