cloudy  trunk
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
save_species.cpp
Go to the documentation of this file.
1 /* This file is part of Cloudy and is copyright (C)1978-2022 by Gary J. Ferland and
2  * others. For conditions of distribution and use see copyright notice in license.txt */
3 /*SaveSpecies generate output for the save species command */
4 #include "cddefines.h"
5 #include "taulines.h"
6 #include "radius.h"
7 #include "save.h"
8 #include "mole.h"
9 #include "mole_priv.h"
10 #include "generic_state.h"
11 #include "iterations.h"
12 #include "dense.h"
13 #include "prt.h"
14 
15 
16 
17 STATIC void SaveSpeciesLines( FILE *ioPUN, const vector<genericState> &speciesList );
18 
19 /* save results for one particular species */
21  const vector<genericState>& SpeciesList,
22  double (*job)(const genericState&),
23  bool lgZonal,
24  const char* chFmt,
25  FILE *ioPUN );
26 
28  const vector<genericState>& SpeciesList,
29  const char *chJob,
30  bool lgZonal,
31  FILE *ioPUN,
32  size_t maxLevels );
33 
34 /*SaveSpecies generate output for the save species command */
36  FILE* ioPUN,
37  long int ipPun)
38 {
39  DEBUG_ENTRY( "SaveSpecies()" );
40 
41  /* first branch; save results for all species if list is empty */
42  vector<genericState> speciesList;
43  if( save.chSaveSpecies[ipPun].size() == 0 )
44  {
45  // loop over species
46  for( size_t i=0; i<mole_global.list.size(); ++i )
47  {
48  speciesList.push_back( genericState(&(mole.species[i])) );
49  }
50  }
51  else
52  {
53  for (vector<string>::iterator name=save.chSaveSpecies[ipPun].begin();
54  name != save.chSaveSpecies[ipPun].end(); ++name)
55  {
56  vector<genericState> v = matchGeneric(name->c_str(),false);
57  if( v.size() == 0 )
58  {
59  fprintf( ioQQQ,"Could not match species '%s', "
60  "use SAVE SPECIES LABELS ALL to get a list of all species."
61  "\nSorry.\n", name->c_str() );
63  }
64  speciesList.insert(speciesList.end(),v.begin(),v.end());
65  }
66  }
67 
68  if( strcmp( save.chSaveArgs[ipPun], "LABE" )==0 )
69  {
70  if( save.lgSaveHeader(ipPun) )
71  {
72  /* save list of species labels */
73  fprintf( ioPUN, "#Species label\tDatabase\n" );
74  for( size_t i=0; i<speciesList.size(); ++i )
75  {
76  fprintf( ioPUN, "%s\t%s\n",
77  speciesList[i].label().c_str(),
78  speciesList[i].database().c_str() );
79  }
80  save.SaveHeaderDone(ipPun);
81  }
82  return;
83  }
84 
85  if( strcmp( save.chSaveArgs[ipPun], "DATA" ) == 0 )
86  {
87  SaveSpeciesLines( ioPUN, speciesList );
88  return;
89  }
90 
91  /* remaining options are column densities, densities, and energies */
92  // max number of levels, for header print
93  size_t mostLevels = 0;
94  for (vector<genericState>::iterator sp=speciesList.begin();
95  sp != speciesList.end(); ++sp)
96  {
97  const molezone *saveSpecies = sp->sp;
98  if( saveSpecies != NULL && saveSpecies != null_molezone &&
99  saveSpecies->levels != NULL )
100  mostLevels = MAX2(mostLevels, saveSpecies->levels->size() );
101  }
102 
103  ASSERT( mostLevels < 10000 );
104 
105  bool lgZonal = true;
106  const char* chJob=NULL, *chFmt=NULL;
107  double (*job)(const genericState&)=NULL;
108  if( strcmp( save.chSaveArgs[ipPun], "ENER" )==0 )
109  {
110  lgZonal = false;
111  chJob = "energies";
112  job = energy;
113  chFmt = "%.5e";
114  }
115  else if( strcmp( save.chSaveArgs[ipPun], "DEPA" )==0 )
116  {
117  lgZonal = true;
118  chJob = "departure coefficients";
119  job = depart;
120  chFmt = "%.5e";
121  }
122  else if( strcmp( save.chSaveArgs[ipPun], "DENS" )==0 )
123  {
124  lgZonal = true;
125  chJob = "densities";
126  chFmt = "%.5e";
127  job = density;
128  }
129  else if( strcmp( save.chSaveArgs[ipPun], "LEVL" )==0 )
130  {
131  lgZonal = true;
132  chJob = "levels";
133  chFmt = "%4.0f";
134  job = levels;
135  }
136  else if( strcmp( save.chSaveArgs[ipPun], "COLU" )==0 )
137  {
138  lgZonal = false;
139  chJob = "column density";
140  chFmt = "%.5e";
141  job = column;
142  }
143  else
144  {
145  fprintf(ioQQQ,"PROBLEM, save species job type \"%s\" not known\n",save.chSaveArgs[ipPun]);
146  TotalInsanity();
147  }
148 
149  if( save.lgSaveHeader(ipPun) )
150  {
151  SaveSpeciesHeader( speciesList, chJob, lgZonal, ioPUN, mostLevels );
152  save.SaveHeaderDone(ipPun);
153  }
154  SaveSpeciesOne( speciesList, job, lgZonal, chFmt, ioPUN );
155 
156  return;
157 }
158 
159 /* print 0.000e+00 as simply 0 */
160 STATIC void PrintShortZero( FILE *ioPUN , const char* chFmt, double arg )
161 {
162  DEBUG_ENTRY( "PrintShortZero()" );
163  if( arg==0. )
164  fprintf(ioPUN,"0");
165  else
166  fprintf(ioPUN,chFmt, arg);
167 
168 }
169 
170 // Switch from initial format to single row per zone
171 const bool lgRowPerZone = true;
172 /* save results for one particular species */
174  const vector<genericState>& speciesList,
175  const char *chJob,
176  bool lgZonal,
177  FILE *ioPUN,
178  size_t maxLevels )
179 {
180  DEBUG_ENTRY( "SaveSpeciesHeader()" );
181 
182  if (!lgRowPerZone)
183  {
184  fprintf( ioPUN, "#%sspecies %s", lgZonal ? "depth\t":"", chJob );
185  for( size_t i = 0; i < maxLevels; ++i )
186  {
187  fprintf( ioPUN, "\t%lu", (unsigned long)i );
188  }
189  fprintf( ioPUN, "\n");
190  }
191  else
192  {
193  int col=0;
194  if (lgZonal)
195  {
196  fprintf( ioPUN, "#depth %s", chJob );
197  ++col;
198  }
199  else
200  {
201  fprintf( ioPUN, "#%s ", chJob );
202  }
203  for (vector<genericState>::const_iterator gs=speciesList.begin();
204  gs != speciesList.end(); ++gs)
205  {
206  if (col == 0)
207  fprintf( ioPUN, "%s",gs->label().c_str() );
208  else
209  fprintf( ioPUN, "\t%s",gs->label().c_str() );
210  ++col;
211  }
212  fprintf( ioPUN,"\n" );
213  }
214 }
215 
217  const vector<genericState>& speciesList,
218  double (*job)(const genericState&),
219  bool lgZonal,
220  const char* chFmt,
221  FILE *ioPUN)
222 {
223  DEBUG_ENTRY( "SaveSpeciesOne()" );
224 
225  int col = 0;
226  if (lgRowPerZone && lgZonal)
227  {
228  fprintf( ioPUN, "%.5e", radius.depth_mid_zone );
229  ++col;
230  }
231  for (vector<genericState>::const_iterator gs=speciesList.begin();
232  gs != speciesList.end(); ++gs)
233  {
234  if (!lgRowPerZone)
235  {
236  if (lgZonal)
237  {
238  fprintf( ioPUN, "%.5e", radius.depth_mid_zone );
239  ++col;
240  }
241  fprintf( ioPUN, "\t%s", gs->label().c_str() );
242  ++col;
243  }
244  if (col > 0)
245  fprintf(ioPUN,"\t");
246  PrintShortZero( ioPUN, chFmt, job(*gs) );
247  ++col;
248  if (!lgRowPerZone)
249  {
250  fprintf(ioPUN,"\n");
251  col = 0;
252  }
253  }
254 
255  if (lgRowPerZone)
256  fprintf(ioPUN,"\n");
257  return;
258 }
259 
261 STATIC void SaveAllSpeciesLabelsLevels( FILE *ioPUN, const vector<genericState> &speciesList )
262 {
263  //Print out number of levels used
264  fprintf( ioPUN, "# The number of levels used in each species.\n" );
265  fprintf( ioPUN, "# Species\tSpectrum\tUsed\tMax.\tDatabase\n" );
266  for( size_t ipSpecies=0; ipSpecies < speciesList.size(); ++ipSpecies )
267  {
268  const molezone *this_mole = speciesList[ ipSpecies ].sp;
269  if( this_mole == NULL || this_mole == null_molezone )
270  continue;
271 
272  species *this_species = (*this_mole).dbase;
273  if( this_species == NULL )
274  continue;
275 
276  fprintf( ioPUN, "%-8s", speciesList[ipSpecies].label().c_str() );
277 
278  string spectralLabel;
279  chemical_to_spectral( speciesList[ ipSpecies ].label(), spectralLabel );
280  fprintf( ioPUN, "\t%s", spectralLabel.c_str() );
281 
282  fprintf( ioPUN, "\t%li", (*this_species).numLevels_local );
283  fprintf( ioPUN, "\t%li", (*this_species).numLevels_max );
284  fprintf( ioPUN, "\t%s", speciesList[ipSpecies].database().c_str() );
285  fprintf( ioPUN, "\n" );
286  }
287 }
288 
289 STATIC void SaveSpeciesLines( FILE *ioPUN, const vector<genericState> &speciesList )
290 {
291  static vector<string> saved_species;
292 
293  vector<genericState> speciesList_local;
294 
295  for( vector<genericState>::const_iterator spcit = speciesList.begin();
296  spcit != speciesList.end(); ++spcit )
297  {
298  bool found = false;
299 
300  for( vector<string>::const_iterator sit = saved_species.begin();
301  sit != saved_species.end(); ++sit )
302  {
303  if( *sit == spcit->label() )
304  {
305  found = true;
306  break;
307  }
308  }
309 
310  if( not found )
311  {
312  speciesList_local.insert( speciesList_local.end(), *spcit );
313  saved_species.push_back( spcit->label() );
314  }
315  }
316 
317  if( speciesList_local.size() == 0 )
318  return;
319 
320  SaveAllSpeciesLabelsLevels( ioPUN, speciesList_local );
321 
322  fprintf( ioPUN,"\n\n");
323 
324  //Species ipLo ipHi gLo gHi wavelen EinA CS Rates
325  fprintf( ioPUN,"Spectrum");
326 
327  if( save.lgSaveDataWn )
328  {
329  fprintf( ioPUN,"\tWavenumbers");
330  }
331  else
332  {
333  fprintf( ioPUN,"\tWavelength");
334  }
335 
336  fprintf( ioPUN,"\tLo");
337  fprintf( ioPUN,"\tHi");
338  fprintf( ioPUN,"\tgLo");
339  fprintf( ioPUN,"\tgHi");
340 
341  if( save.lgSaveDataGf )
342  {
343  fprintf( ioPUN,"\t gf ");
344  }
345  else
346  {
347  fprintf( ioPUN,"\tEinstein A");
348  }
349  fprintf( ioPUN,"\tColl_Str");
350  fprintf( ioPUN,"\tgbar");
351 
352  if( save.lgSaveDataRates )
353  {
354  fprintf( ioPUN,"\tRate electron");
355  fprintf( ioPUN,"\tRate proton");
356  fprintf( ioPUN,"\tRate He+ ");
357  fprintf( ioPUN,"\tRate Alpha");
358  fprintf( ioPUN,"\tRate Atom H");
359  fprintf( ioPUN,"\tRate Atom He");
360  fprintf( ioPUN,"\tRate Ortho");
361  fprintf( ioPUN,"\tRate Para");
362  fprintf( ioPUN,"\tRate H2");
363  }
364 
365  fprintf( ioPUN,"\n");
366 
367 
368  for( size_t ipSpecies=0; ipSpecies < speciesList_local.size(); ++ipSpecies )
369  {
370  const molezone *this_mole = speciesList_local[ ipSpecies ].sp;
371  if( this_mole == NULL || this_mole == null_molezone )
372  continue;
373  if( (*this_mole).lines == NULL )
374  continue;
375 
376  species *this_species = (*this_mole).dbase;
377  if( this_species == NULL )
378  continue;
379 
380  for (TransitionList::iterator tr = (*this_mole).lines->begin();
381  tr != (*this_mole).lines->end(); ++tr )
382  {
383  long ipLo = tr->ipLo() +1;
384  long ipHi = tr->ipHi() +1;
385  int nelem = tr->Hi()->nelem();
386 
387  if( nelem != -1 && !dense.lgElmtOn[nelem-1] )
388  {
389  continue;
390  }
391 
392  if( ipHi >= (*this_species).numLevels_local )
393  {
394  continue;
395  }
396 
397  string spectralLabel;
398  chemical_to_spectral( speciesList_local[ ipSpecies ].label(), spectralLabel );
399  fprintf( ioPUN,"%-8s", spectralLabel.c_str() );
400 
401  if( save.lgSaveDataWn )
402  {
403  fprintf( ioPUN,"\t%.3e", tr->EnergyWN() );
404  }
405  else
406  {
407  fprintf( ioPUN, "\t" );
408  prt_wl( ioPUN, realnum( tr->WLAng() ) );
409  }
410 
411  fprintf( ioPUN,"\t%li", ipLo);
412  fprintf( ioPUN,"\t%li", ipHi);
413  fprintf( ioPUN,"\t%i", int( tr->Lo()->g() ) );
414  fprintf( ioPUN,"\t%i", int( tr->Hi()->g() ) );
415 
416  if( save.lgSaveDataGf )
417  {
418  fprintf( ioPUN,"\t%.3e", tr->Emis().gf() );
419  }
420  else
421  {
422  fprintf( ioPUN,"\t%.3e", tr->Emis().Aul() );
423  }
424 
425  fprintf( ioPUN,"\t%.3e", tr->Coll().col_str());
426  fprintf( ioPUN,"\t%i", tr->Coll().is_gbar() );
427 
428  if( save.lgSaveDataRates )
429  {
430  for( long intCollNo=0; intCollNo<ipNCOLLIDER; intCollNo++)
431  {
432  fprintf( ioPUN,"\t%.3e",tr->Coll().rate_coef_ul()[intCollNo]);
433  }
434  }
435 
436  fprintf( ioPUN,"\n");
437  }
438  }
439 }
440 
441 void SaveSpeciesOptDep( const long int ipPun, const string &speciesLabel )
442 {
443  DEBUG_ENTRY( "SaveSpeciesOptDep()" );
444 
445  if( ! iterations.lgLastIt )
446  return;
447 
448  if( save.lgSaveHeader(ipPun) )
449  {
450  fprintf( save.params[ipPun].ipPnunit,
451  "#hi\tlo\tWavelength(A)\ttau\n" );
452  save.SaveHeaderDone( ipPun );
453  }
454 
456  getSpecies( speciesLabel, species );
457 
458  if( species.sp->lines == NULL )
459  {
460  fprintf( ioQQQ,
461  "WARNING: Species '%s' does not have any data for 'save species optical depth'.\n",
462  species.label().c_str() );
463  return;
464  }
465 
466  for( TransitionList::iterator tr = (*species.sp->lines).begin();
467  tr != (*species.sp->lines).end(); ++tr)
468  {
469  if( (*tr).ipCont() <= 0 )
470  continue;
471 
472  fprintf( save.params[ipPun].ipPnunit,
473  "%i\t%i\t%.5e\t%.5e\n",
474  (*tr).ipHi()+1,
475  (*tr).ipLo()+1,
476  (*tr).WLAng(),
477  (*tr).Emis().TauIn() * SQRTPI );
478  }
479 }
480 
481 class Field
482 {
483 public:
484  const char* label;
485  const char* fmt;
486 private:
487  double m_value;
488 public:
489  Field(const char* label, const char* fmt, double value) :
490  label(label), fmt(fmt), m_value(value)
491  {}
492  double value() const
493  {
494  return m_value;
495  }
496 };
497 
498 STATIC void doHeader(FILE *punit, const Field& f)
499 {
500  DEBUG_ENTRY( "doHeader()" );
501  fprintf(punit,"%s",f.label);
502 }
503 STATIC void doData(FILE *punit, const Field& f)
504 {
505  DEBUG_ENTRY( "doData()" );
506  fprintf(punit,f.fmt,f.value());
507 }
508 
509 STATIC inline bool isCatalystReactant(const mole_reaction& rate, int i)
510 {
511  return rate.rvector[i]!=NULL;
512 }
513 STATIC inline bool isCatalystProduct(const mole_reaction& rate, int i)
514 {
515  return rate.pvector[i]!=NULL;
516 }
517 STATIC inline bool isDestroyed(const mole_reaction& rate, int i)
518 {
519  return !isCatalystReactant(rate,i) && rate.rvector_excit[i]==NULL;
520 }
521 STATIC inline bool isCreated(const mole_reaction& rate, int i)
522 {
523  return !isCatalystProduct(rate,i) && rate.pvector_excit[i]==NULL;
524 }
525 
526 /* Save all rates involving specified species */
527 void mole_save(FILE *punit, const char speciesname[], const char args[], bool lgHeader, bool lgData, bool lgCoef, double depth)
528 {
529  DEBUG_ENTRY( "mole_save()" );
530 
531  molecule *sp = findspecies(speciesname);
532 
533  if( sp == null_mole )
534  {
535  fprintf( punit, " Species %s could not be found. Note that labels are case-sensitive in this context.\n", speciesname );
536  return;
537  }
538 
539  void (*doTask)(FILE *punit, const Field& f);
540 
541  if (lgHeader)
542  doTask=doHeader;
543  else
544  doTask=doData;
545 
546  doTask(punit,Field("Depth","%.5e",depth));
547 
548  for (mole_reaction_i p
549  =mole_priv::reactab.begin(); p != mole_priv::reactab.end(); ++p)
550  {
551  const mole_reaction &rate = *p->second;
552  bool ipthis = false;
553 
554  for (int i=0;i<rate.nreactants && !ipthis;i++)
555  {
556  if( rate.reactants[i] == sp )
557  {
558  if( ( strcmp( args, "DEST" )==0 && isDestroyed(rate,i) ) ||
559  ( strcmp( args, "DFLT" )==0 && isDestroyed(rate,i) ) ||
560  ( strcmp( args, "CATA" )==0 && isCatalystReactant(rate,i) ) ||
561  strcmp( args, "ALL " )==0 )
562  ipthis = true;
563  }
564  }
565 
566  for(int i=0;i<rate.nproducts && !ipthis;i++)
567  {
568  if( rate.products[i] == sp )
569  {
570  if( ( strcmp( args, "CREA" )==0 && isCreated(rate,i) ) ||
571  ( strcmp( args, "DFLT" )==0 && isCreated(rate,i) ) ||
572  ( strcmp( args, "CATA" )==0 && isCatalystProduct(rate,i) ) ||
573  strcmp( args, "ALL " )==0 )
574  ipthis = true;
575  }
576  }
577 
578  if(ipthis)
579  {
580  double ratevi=0.0;
581  if (lgData)
582  {
583  ratevi = mole.reaction_rks[ rate.index ];
584  if( !lgCoef )
585  {
586  for(int i=0;i<rate.nreactants;i++)
587  {
588  ratevi *= mole.species[ rate.reactants[i]->index ].den;
589  }
590  }
591  }
592  fprintf(punit,"\t");
593 
594  doTask(punit,Field(rate.label.c_str(),"%.3e",ratevi));
595  }
596  }
597  fprintf(punit,"\n");
598 }
599 
600 // Helper class for sorting rates
601 class RateCmp
602 {
603  const vector<double>& m_stack;
604 public:
605  RateCmp(const vector<double>& stack) : m_stack(stack) {}
606  bool operator()(size_t a, size_t b)
607  {
608  return m_stack[a] > m_stack[b];
609  }
610 };
611 
612 void mole_dominant_rates( const vector<const molecule*>& debug_list,
613  FILE *ioOut,
614  bool lgPrintReagents, size_t NPRINT, double fprint )
615 {
616  vector<double> snkx, srcx;
617  vector<mole_reaction *> ratesnk, ratesrc;
618 
619  double src_total = 0.0, snk_total = 0.0;
620 
621  for(mole_reaction_i p=mole_priv::reactab.begin();
622  p != mole_priv::reactab.end(); ++p)
623  {
624  mole_reaction *rate = &(*p->second);
625  double rk = mole.reaction_rks[ rate->index ];
626 
627  double rate_tot = rk;
628  for(long i=0;i<rate->nreactants;i++)
629  {
630  rate_tot *= mole.species[ rate->reactants[i]->index ].den;
631  }
632 
633  int nrate=0;
634  for(size_t s=0; s<debug_list.size(); ++s)
635  {
636  for(long i=0;i<rate->nproducts;++i)
637  {
638  molecule *sp = rate->products[i];
639  if (sp == debug_list[s] && rate->pvector[i] == NULL)
640  {
641  ++nrate;
642  }
643  }
644  for(long i=0;i<rate->nreactants;++i)
645  {
646  molecule *sp = rate->reactants[i];
647  if (sp == debug_list[s] && rate->rvector[i] == NULL)
648  {
649  --nrate;
650  }
651  }
652  }
653  if (nrate > 0)
654  {
655  srcx.push_back(nrate*rate_tot);
656  ratesrc.push_back(rate);
657  src_total += nrate*rate_tot;
658  }
659  else if (nrate < 0)
660  {
661  snkx.push_back(-nrate*rate_tot);
662  ratesnk.push_back(rate);
663  snk_total -= nrate*rate_tot;
664  }
665  }
666 
667  if (!ratesrc.empty())
668  {
669  vector<size_t> isrc(ratesrc.size());
670  for (size_t i=0; i<ratesrc.size(); ++i)
671  isrc[i] = i;
672  sort(isrc.begin(),isrc.end(),RateCmp(srcx));
673 
674  fprintf( ioOut, "Src %13.7g: ",src_total);
675  for (size_t i=0; i<ratesrc.size(); ++i)
676  {
677  if (i == NPRINT || srcx[isrc[i]] < fprint*src_total)
678  {
679  break;
680  }
681  fprintf( ioOut, "%20.20s %13.7g",
682  ratesrc[isrc[i]]->label.c_str(),srcx[isrc[i]]);
683  if (lgPrintReagents)
684  {
685  fprintf( ioOut, " [");
686  for (long j=0;j<ratesrc[isrc[i]]->nreactants;j++)
687  {
688  if (j)
689  {
690  fprintf( ioOut, "," );
691  }
692  fprintf( ioOut, "%-6.6s %13.7g",
693  ratesrc[isrc[i]]->reactants[j]->label.c_str(),
694  mole.species[ ratesrc[isrc[i]]->reactants[j]->index ].den);
695  }
696  fprintf( ioOut, "]" );
697  }
698  else
699  {
700  fprintf( ioOut, ";" );
701  }
702  }
703  }
704  if (!ratesnk.empty())
705  {
706  vector<size_t> isnk(ratesnk.size());
707  for (size_t i=0; i<ratesnk.size(); ++i)
708  isnk[i] = i;
709  sort(isnk.begin(),isnk.end(),RateCmp(snkx));
710 
711  fprintf( ioOut, " Snk %13.7g: ", snk_total);
712  for (size_t i=0; i<ratesnk.size(); ++i)
713  {
714  if (i == NPRINT || snkx[isnk[i]] < fprint*snk_total)
715  {
716  break;
717  }
718  fprintf( ioOut, "%20.20s %13.7g",
719  ratesnk[isnk[i]]->label.c_str(), snkx[isnk[i]] );
720  if (lgPrintReagents)
721  {
722  fprintf( ioOut, " [" );
723  for (long j=0;j<ratesnk[isnk[i]]->nreactants;j++)
724  {
725  if (j)
726  {
727  fprintf( ioOut, "," );
728  }
729  fprintf( ioOut, "%-6.6s %13.7g",
730  ratesnk[isnk[i]]->reactants[j]->label.c_str(),
731  mole.species[ ratesnk[isnk[i]]->reactants[j]->index ].den);
732  }
733  fprintf( ioOut, "]" );
734  }
735  else
736  {
737  fprintf( ioOut, ";" );
738  }
739  }
740  }
741  fprintf( ioOut, "\n" );
742 
743  return;
744 }
745 
746 void mole_print_species_reactions( molecule *speciesToPrint )
747 {
748  if( speciesToPrint==NULL )
749  {
750  fprintf( ioQQQ,"\n NULL species found in mole_print_species_reactions.\n" );
751  return;
752  }
753 
754  long numReacts = 0;
755 
756  fprintf( ioQQQ,"\n Reactions involving species %s:\n", speciesToPrint->label.c_str() );
757  for( mole_reaction_i p=mole_priv::reactab.begin(); p!=mole_priv::reactab.end(); ++p )
758  {
759  mole_reaction &rate = *p->second;
760  for( long i=0; i<rate.nreactants; i++ )
761  {
762  molecule *sp = rate.reactants[i];
763  if(rate.rvector[i]==NULL && sp==speciesToPrint )
764  {
765  double drate = mole.reaction_rks[ rate.index ];
766  for (long j=0; j<rate.nreactants; ++j)
767  {
768  drate *= mole.species[ rate.reactants[j]->index ].den;
769  }
770  fprintf( ioQQQ,"%s rate = %g\n", rate.label.c_str(), drate );
771  numReacts++;
772  }
773  }
774  for( long i=0; i<rate.nproducts; i++ )
775  {
776  molecule *sp = rate.products[i];
777  if(rate.pvector[i]==NULL && sp==speciesToPrint )
778  {
779  double drate = mole.reaction_rks[ rate.index ];
780  for (long j=0; j<rate.nreactants; ++j)
781  {
782  drate *= mole.species[ rate.reactants[j]->index ].den;
783  }
784  fprintf( ioQQQ,"%s rate = %g\n", rate.label.c_str(), drate );
785  numReacts++;
786  }
787  }
788  }
789  fprintf( ioQQQ," End of reactions involving species %s. There were %li.\n", speciesToPrint->label.c_str(), numReacts );
790 
791  return;
792 }
molecule * reactants[MAXREACTANTS]
Definition: mole_priv.h:53
t_mole_global mole_global
Definition: mole.cpp:7
vector< double > reaction_rks
Definition: mole.h:470
void prt_wl(FILE *ioOUT, realnum wl)
Definition: prt.cpp:44
STATIC void PrintShortZero(FILE *ioPUN, const char *chFmt, double arg)
STATIC long int ipPun
Definition: save_do.cpp:367
molecule * null_mole
NORETURN void TotalInsanity(void)
Definition: service.cpp:971
double depart(const genericState &gs)
STATIC void doHeader(FILE *punit, const Field &f)
STATIC bool isCreated(const mole_reaction &rate, int i)
double value() const
int nreactants
Definition: mole_priv.h:52
const vector< double > & m_stack
map< string, count_ptr< mole_reaction > > reactab
vector< genericState > matchGeneric(const char *chLabel, bool lgValidate)
void getSpecies(const string &speciesLabel, genericState &species)
TransitionList * lines
Definition: mole.h:418
STATIC bool isCatalystProduct(const mole_reaction &rate, int i)
STATIC void SaveSpeciesHeader(const vector< genericState > &SpeciesList, const char *chJob, bool lgZonal, FILE *ioPUN, size_t maxLevels)
FILE * ioQQQ
Definition: cddefines.cpp:7
FILE * ipPnunit
Definition: save.h:195
Definition: mole.h:378
Definition: mole.h:142
STATIC void SaveSpeciesOne(const vector< genericState > &SpeciesList, double(*job)(const genericState &), bool lgZonal, const char *chFmt, FILE *ioPUN)
double m_value
const char * label
void mole_save(FILE *punit, const char speciesname[], const char args[], bool lgHeader, bool lgData, bool lgCoef, double depth)
t_dense dense
Definition: global.cpp:15
void SaveAllSpeciesLabelsLevels(FILE *ioPUN)
RateCmp(const vector< double > &stack)
void SaveSpeciesOptDep(const long int ipPun, const string &speciesLabel)
molecule * products[MAXPRODUCTS]
Definition: mole_priv.h:56
const char * fmt
qList * levels
Definition: mole.h:417
string label
Definition: mole_priv.h:51
STATIC bool isCatalystReactant(const mole_reaction &rate, int i)
molecule * rvector[MAXREACTANTS]
Definition: mole_priv.h:54
double energy(const genericState &gs)
#define STATIC
Definition: cddefines.h:118
char chSaveArgs[LIMPUN][5]
Definition: save.h:383
void chemical_to_spectral(const string chLabelChem, string &chLabelSpec)
Definition: species.cpp:983
t_mole_local mole
Definition: mole.cpp:8
map< string, count_ptr< mole_reaction > >::iterator mole_reaction_i
Definition: mole_priv.h:38
molecule * findspecies(const char buf[])
float realnum
Definition: cddefines.h:124
valarray< class molezone > species
Definition: mole.h:468
#define EXIT_FAILURE
Definition: cddefines.h:168
STATIC void SaveSpeciesLines(FILE *ioPUN, const vector< genericState > &speciesList)
bool lgSaveDataRates
Definition: save.h:506
size_t size() const
Definition: quantumstate.h:121
#define cdEXIT(FAIL)
Definition: cddefines.h:482
int index
Definition: mole.h:194
Field(const char *label, const char *fmt, double value)
double depth_mid_zone
Definition: radius.h:31
t_iterations iterations
Definition: iterations.cpp:6
double column(const genericState &gs)
const molezone * sp
Definition: generic_state.h:19
t_radius radius
Definition: radius.cpp:5
SaveParams params[LIMPUN]
Definition: save.h:278
bool lgElmtOn[LIMELM]
Definition: dense.h:160
STATIC void doData(FILE *punit, const Field &f)
double levels(const genericState &gs)
void SaveHeaderDone(int ipPun)
Definition: save.h:364
#define ASSERT(exp)
Definition: cddefines.h:613
bool lgLastIt
Definition: iterations.h:47
molecule * pvector[MAXPRODUCTS]
Definition: mole_priv.h:57
void mole_dominant_rates(const vector< const molecule * > &debug_list, FILE *ioOut, bool lgPrintReagents, size_t NPRINT, double fprint)
double density(const genericState &gs)
bool lgSaveDataGf
Definition: save.h:506
#define DEBUG_ENTRY(funcname)
Definition: cddefines.h:723
const bool lgRowPerZone
string label
Definition: mole.h:156
void SaveSpecies(FILE *ioPUN, long int ipPun)
#define MAX2(a, b)
Definition: cddefines.h:824
int fprintf(const Output &stream, const char *format,...)
Definition: service.cpp:1121
MoleculeList list
Definition: mole.h:365
vector< string > chSaveSpecies[LIMPUN]
Definition: save.h:385
molecule * rvector_excit[MAXREACTANTS]
Definition: mole_priv.h:55
string label() const
bool lgSaveDataWn
Definition: save.h:506
t_save save
Definition: save.cpp:5
STATIC bool isDestroyed(const mole_reaction &rate, int i)
bool lgSaveHeader(int ipPun) const
Definition: save.h:360
molezone * null_molezone
bool operator()(size_t a, size_t b)
void mole_print_species_reactions(molecule *speciesToPrint)
molecule * pvector_excit[MAXPRODUCTS]
Definition: mole_priv.h:58