00001
00002
00003
00004
00005 #include "cddefines.h"
00006 #include "physconst.h"
00007 #include "mole.h"
00008 #include "mole_co_priv.h"
00009 #include "hmi.h"
00010 #include "rfield.h"
00011 #include "dense.h"
00012 #include "ionbal.h"
00013 #include "grainvar.h"
00014 #include "timesc.h"
00015
00016
00017
00018
00019
00020 enum spectype {MOLECULE, OTHER};
00021 enum molstate {ACTIVE, PASSIVE};
00022
00023 STATIC void newelement(const char label[], int ipion,
00024 int priority);
00025 STATIC struct molecule *newspecies(const char label[7], enum spectype type,
00026 enum molstate state, realnum *location, double frac0);
00027 STATIC struct chem_element_s *findelement(const char buf[]);
00028 STATIC int isactive(data_u *dat);
00029 STATIC int ispassive(data_u *dat);
00030 STATIC int isCOnet(data_u *dat);
00031
00032 struct chem_element_s **chem_element;
00033
00034
00035 struct chem_element_s *element_list[LIMELM];
00036 int32 *ipiv;
00037 realnum *tot_ion;
00038
00039 struct mole_priv_s mole_priv;
00040 struct molecule null_mole;
00041
00042
00043
00044 void CO_Init(void)
00045 {
00046 long int i,
00047 nelem;
00048 struct molecule *sp;
00049 static bool lgCO_Init_called=false;
00050
00051 DEBUG_ENTRY( "CO_Init()" );
00052
00053
00054
00055 co.co_nzone = -2;
00056 co.iteration_co = -2;
00057
00058
00059
00060
00061 if( lgCO_Init_called )
00062 {
00063 return;
00064 }
00065
00066
00067 lgCO_Init_called = true;
00068
00069 mole_priv.spectab = newhash(NULL);
00070 mole_priv.reactab = newhash(NULL);
00071 mole_priv.elemtab = newhash(NULL);
00072
00073 for(nelem=ipHYDROGEN; nelem<LIMELM; ++nelem )
00074 {
00075 element_list[nelem] = NULL;
00076 mole.lgElem_in_chemistry[nelem] = false;
00077 }
00078
00079
00080
00081
00082
00083
00084
00085 newelement("C" ,ipCARBON,3);
00086 newelement("^13C",ipCARBON,0);
00087 newelement("O" ,ipOXYGEN,2);
00088 newelement("Si",ipSILICON,5);
00089 newelement("N" ,ipNITROGEN,4);
00090 newelement("S" ,ipSULPHUR,6);
00091 newelement("Cl",ipCHLORINE,7);
00092 newelement("H" ,ipHYDROGEN,1);
00093 newelement("He",ipHELIUM,0);
00094 newelement("Mg",ipMAGNESIUM,0);
00095 newelement("Fe",ipIRON,0);
00096
00097
00098
00099
00100
00101
00102
00103
00104
00105
00106 mole.num_comole_calc = mole.num_comole_tot = mole.num_elements = 0;
00107 null_mole.hevmol = null_mole.hevcol = 0.;
00108 null_mole.index = -1;
00109
00110 sp = newspecies("C ",MOLECULE,ACTIVE,NULL,1.00e+00);
00111 sp = newspecies("O ",MOLECULE,ACTIVE,NULL,1.13e+05);
00112 sp = newspecies("Si ",MOLECULE,ACTIVE,NULL,3.56e-04);
00113 sp = newspecies("N ",MOLECULE,ACTIVE,NULL,2.25e+00);
00114 sp = newspecies("S ",MOLECULE,ACTIVE,NULL,6.90e-03);
00115 sp = newspecies("Cl ",MOLECULE,ACTIVE,NULL,6.90e-03);
00116 sp = newspecies("C+ ",MOLECULE,ACTIVE,NULL,6.79e+04);
00117 sp = newspecies("O+ ",MOLECULE,ACTIVE,NULL,1.58e+01);
00118 sp = newspecies("Si+ ",MOLECULE,ACTIVE,NULL,1.79e+02);
00119 sp = newspecies("N+ ",MOLECULE,ACTIVE,NULL,2.62e-10);
00120 sp = newspecies("S+ ",MOLECULE,ACTIVE,NULL,1.79e+03);
00121 sp = newspecies("Cl+ ",MOLECULE,ACTIVE,NULL,3.71e-09);
00122 sp = newspecies("CH ",MOLECULE,ACTIVE,NULL,1.10e-06);
00123 sp = newspecies("CH+ ",MOLECULE,ACTIVE,NULL,6.99e-05);
00124 sp = newspecies("OH ",MOLECULE,ACTIVE,NULL,2.15e-03);
00125 sp = newspecies("OH+ ",MOLECULE,ACTIVE,NULL,2.45e-04);
00126 sp = newspecies("O2 ",MOLECULE,ACTIVE,NULL,1.91e-07);
00127 sp = newspecies("CO ",MOLECULE,ACTIVE,NULL,4.05e-05);
00128 sp = newspecies("^13CO ",MOLECULE,ACTIVE,NULL,4.05e-05);
00129 sp = newspecies("CO+ ",MOLECULE,ACTIVE,NULL,2.16e-06);
00130 sp = newspecies("H2O ",MOLECULE,ACTIVE,NULL,1.60e-11);
00131 sp = newspecies("H2O+ ",MOLECULE,ACTIVE,NULL,1.27e-10);
00132 sp = newspecies("O2+ ",MOLECULE,ACTIVE,NULL,1.17e-06);
00133 sp = newspecies("H3O+ ",MOLECULE,ACTIVE,NULL,3.44e-17);
00134 sp = newspecies("CH2+ ",MOLECULE,ACTIVE,NULL,3.71e-09);
00135 sp = newspecies("CH2 ",MOLECULE,ACTIVE,NULL,8.59e-13);
00136 sp = newspecies("HCO+ ",MOLECULE,ACTIVE,NULL,4.01e-10);
00137 sp = newspecies("CH3+ ",MOLECULE,ACTIVE,NULL,7.02e-16);
00138 sp = newspecies("CH3 ",MOLECULE,ACTIVE,NULL,4.59e-19);
00139 sp = newspecies("CH4 ",MOLECULE,ACTIVE,NULL,9.12e-28);
00140 sp = newspecies("CH4+ ",MOLECULE,ACTIVE,NULL,1.03e-28);
00141 sp = newspecies("CH5+ ",MOLECULE,ACTIVE,NULL,8.16e-28);
00142 sp = newspecies("SiH2+ ",MOLECULE,ACTIVE,NULL,2.07e-13);
00143 sp = newspecies("SiH ",MOLECULE,ACTIVE,NULL,1.80e-14);
00144 sp = newspecies("SiOH+ ",MOLECULE,ACTIVE,NULL,5.39e-15);
00145 sp = newspecies("SiO ",MOLECULE,ACTIVE,NULL,3.89e-14);
00146 sp = newspecies("SiO+ ",MOLECULE,ACTIVE,NULL,2.27e-08);
00147 sp = newspecies("N2 ",MOLECULE,ACTIVE,NULL,2.46e-17);
00148 sp = newspecies("N2+ ",MOLECULE,ACTIVE,NULL,2.22e-17);
00149 sp = newspecies("NO ",MOLECULE,ACTIVE,NULL,5.99e-12);
00150 sp = newspecies("NO+ ",MOLECULE,ACTIVE,NULL,1.31e-11);
00151 sp = newspecies("S2 ",MOLECULE,ACTIVE,NULL,1.21e-11);
00152 sp = newspecies("S2+ ",MOLECULE,ACTIVE,NULL,1.00e-13);
00153 sp = newspecies("OCN ",MOLECULE,ACTIVE,NULL,3.86e-11);
00154 sp = newspecies("OCN+ ",MOLECULE,ACTIVE,NULL,1.65e-22);
00155 sp = newspecies("NH ",MOLECULE,ACTIVE,NULL,1.30e-09);
00156 sp = newspecies("NH+ ",MOLECULE,ACTIVE,NULL,2.20e-10);
00157 sp = newspecies("NH2 ",MOLECULE,ACTIVE,NULL,6.78e-20);
00158 sp = newspecies("NH2+ ",MOLECULE,ACTIVE,NULL,1.71e-16);
00159 sp = newspecies("NH3 ",MOLECULE,ACTIVE,NULL,1.25e-29);
00160 sp = newspecies("NH3+ ",MOLECULE,ACTIVE,NULL,3.03e-23);
00161 sp = newspecies("NH4+ ",MOLECULE,ACTIVE,NULL,1.15e-33);
00162 sp = newspecies("CN ",MOLECULE,ACTIVE,NULL,3.38e-11);
00163 sp = newspecies("CN+ ",MOLECULE,ACTIVE,NULL,5.07e-11);
00164 sp = newspecies("HCN ",MOLECULE,ACTIVE,NULL,1.07e-17);
00165 sp = newspecies("HCN+ ",MOLECULE,ACTIVE,NULL,9.89e-17);
00166 sp = newspecies("HNO ",MOLECULE,ACTIVE,NULL,9.09e-21);
00167 sp = newspecies("HNO+ ",MOLECULE,ACTIVE,NULL,1.91e-18);
00168 sp = newspecies("HS ",MOLECULE,ACTIVE,NULL,1.71e-14);
00169 sp = newspecies("HS+ ",MOLECULE,ACTIVE,NULL,4.83e-13);
00170 sp = newspecies("CS ",MOLECULE,ACTIVE,NULL,6.69e-18);
00171 sp = newspecies("CS+ ",MOLECULE,ACTIVE,NULL,4.12e-11);
00172 sp = newspecies("NO2 ",MOLECULE,ACTIVE,NULL,2.67e-26);
00173 sp = newspecies("NO2+ ",MOLECULE,ACTIVE,NULL,2.41e-23);
00174 sp = newspecies("NS ",MOLECULE,ACTIVE,NULL,3.12e-11);
00175 sp = newspecies("NS+ ",MOLECULE,ACTIVE,NULL,6.81e-13);
00176 sp = newspecies("SO ",MOLECULE,ACTIVE,NULL,4.00e-15);
00177 sp = newspecies("SO+ ",MOLECULE,ACTIVE,NULL,1.20e-07);
00178 sp = newspecies("SiN ",MOLECULE,ACTIVE,NULL,7.03e-12);
00179 sp = newspecies("SiN+ ",MOLECULE,ACTIVE,NULL,7.60e-14);
00180 sp = newspecies("N2O ",MOLECULE,ACTIVE,NULL,1.05e-11);
00181 sp = newspecies("HCS+ ",MOLECULE,ACTIVE,NULL,8.91e-17);
00182 sp = newspecies("OCS ",MOLECULE,ACTIVE,NULL,8.83e-24);
00183 sp = newspecies("OCS+ ",MOLECULE,ACTIVE,NULL,1.72e-21);
00184 sp = newspecies("HNC ",MOLECULE,ACTIVE,NULL,4.00e-15);
00185 sp = newspecies("HCNH+ ",MOLECULE,ACTIVE,NULL,4.00e-15);
00186 sp = newspecies("C2 ",MOLECULE,ACTIVE,NULL,3.71e-09);
00187 sp = newspecies("C2+ ",MOLECULE,ACTIVE,NULL,8.59e-13);
00188 sp = newspecies("CCl ",MOLECULE,ACTIVE,NULL,4.00e-15);
00189 sp = newspecies("ClO ",MOLECULE,ACTIVE,NULL,4.00e-15);
00190 sp = newspecies("HCl+ ",MOLECULE,ACTIVE,NULL,4.00e-15);
00191 sp = newspecies("HCl ",MOLECULE,ACTIVE,NULL,4.00e-15);
00192 sp = newspecies("H2Cl+ ",MOLECULE,ACTIVE,NULL,4.00e-15);
00193 sp = newspecies("CCl+ ",MOLECULE,ACTIVE,NULL,4.00e-15);
00194 sp = newspecies("H2CCl+",MOLECULE,ACTIVE,NULL,0.00e+00);
00195 sp = newspecies("ClO+ ",MOLECULE,ACTIVE,NULL,4.00e-15);
00196
00197 if(gv.lgDustOn() && mole.lgGrain_mole_deplete )
00198 {
00199 sp = newspecies("COgrn ",MOLECULE,ACTIVE,NULL,1.00e-15);
00200 sp = newspecies("H2Ogrn",MOLECULE,ACTIVE,NULL,0.00e+00);
00201 sp = newspecies("OHgrn ",MOLECULE,ACTIVE,NULL,1.00e-15);
00202 }
00203 sp = newspecies("C2H ",MOLECULE,ACTIVE,NULL,4.00e-15);
00204 sp = newspecies("C2H+ ",MOLECULE,ACTIVE,NULL,4.00e-15);
00205 sp = newspecies("C3 ",MOLECULE,ACTIVE,NULL,1.00e-15);
00206 sp = newspecies("C3+ ",MOLECULE,ACTIVE,NULL,4.00e-15);
00207 sp = newspecies("C3H+ ",MOLECULE,ACTIVE,NULL,4.00e-15);
00208 sp = newspecies("C3H ",MOLECULE,ACTIVE,NULL,4.00e-15);
00209 sp = newspecies("C2H2 ",MOLECULE,ACTIVE,NULL,4.00e-15);
00210 sp = newspecies("C2H2+ ",MOLECULE,ACTIVE,NULL,4.00e-15);
00211 sp = newspecies("C2H3+ ",MOLECULE,ACTIVE,NULL,4.00e-15);
00212 sp = newspecies("N2H+ ",MOLECULE,ACTIVE,NULL,4.00e-15);
00213
00214 sp = newspecies("e- ",OTHER,PASSIVE,&(dense.eden_f),0.00e+00);
00215 sp->nElec = -1; sp->mole_mass = (realnum)ELECTRON_MASS;
00216
00217 sp = newspecies("Fe ",MOLECULE,PASSIVE,&(dense.xIonDense[ipIRON][0]),0.00e+00);
00218 sp = newspecies("Fe+ ",MOLECULE,PASSIVE,&(dense.xIonDense[ipIRON][1]),0.00e+00);
00219 sp = newspecies("H ",MOLECULE,PASSIVE,&(dense.xIonDense[ipHYDROGEN][0]),0.00e+00);
00220 sp = newspecies("H- ",MOLECULE,PASSIVE,&(hmi.Hmolec[ipMHm]),0.00e+00);
00221 sp = newspecies("H+ ",MOLECULE,PASSIVE,&(dense.xIonDense[ipHYDROGEN][1]),0.00e+00);
00222 sp = newspecies("H2 ",MOLECULE,PASSIVE,&(hmi.Hmolec[ipMH2g]),0.00e+00);
00223 sp = newspecies("H2* ",MOLECULE,PASSIVE,&(hmi.Hmolec[ipMH2s]),0.00e+00);
00224 sp = newspecies("H2+ ",MOLECULE,PASSIVE,&(hmi.Hmolec[ipMH2p]),0.00e+00);
00225 sp = newspecies("H3+ ",MOLECULE,PASSIVE,&(hmi.Hmolec[ipMH3p]),0.00e+00);
00226 sp = newspecies("He ",MOLECULE,PASSIVE,&(dense.xIonDense[ipHELIUM][0]),0.00e+00);
00227 sp = newspecies("He+ ",MOLECULE,PASSIVE,&(dense.xIonDense[ipHELIUM][1]),0.00e+00);
00228 sp = newspecies("Mg ",MOLECULE,PASSIVE,&(dense.xIonDense[ipMAGNESIUM][0]),0.00e+00);
00229 sp = newspecies("Mg+ ",MOLECULE,PASSIVE,&(dense.xIonDense[ipMAGNESIUM][1]),0.00e+00);
00230
00231
00232 COmole = (struct molecule **)MALLOC((size_t)mole.num_comole_tot*
00233 sizeof(struct molecule *));
00234
00235
00236 i = makeplist(mole_priv.spectab,(void **)COmole,
00237 mole.num_comole_calc,isactive);
00238 ASSERT (i == mole.num_comole_calc);
00239
00240
00241 i = makeplist(mole_priv.spectab,(void **)COmole+mole.num_comole_calc,
00242 mole.num_comole_tot-mole.num_comole_calc,ispassive);
00243 ASSERT (i == mole.num_comole_tot-mole.num_comole_calc);
00244
00245
00246 for(i=0;i<mole.num_comole_tot;i++)
00247 {
00248 COmole[i]->index = i;
00249 }
00250
00251
00252 for(i=0;i<mole.num_comole_calc;i++)
00253 {
00254 sp = COmole[i];
00255 if(sp->active && sp->n_nuclei == 1)
00256 {
00257 if(sp->nElec == 0)
00258 {
00259 element_list[sp->nelem_hevmol]->ipMl = i;
00260 }
00261 else if(sp->nElec == 1)
00262 {
00263 element_list[sp->nelem_hevmol]->ipMlP = i;
00264 }
00265 }
00266 }
00267
00268
00269 chem_element = (struct chem_element_s **)
00270 MALLOC((size_t)(mole.num_elements)*sizeof(struct chem_element_s *));
00271 i = makeplist(mole_priv.elemtab,(void **)chem_element,mole.num_elements,isCOnet);
00272 ASSERT(i == mole.num_elements);
00273
00274
00275
00276 mole.amat = (double **)MALLOC((size_t)mole.num_comole_calc*sizeof(double *));
00277 mole.amat[0] = (double *)MALLOC((size_t)mole.num_comole_calc*
00278 mole.num_comole_calc*sizeof(double));
00279 for(i=1;i<mole.num_comole_calc;i++)
00280 {
00281 mole.amat[i] = mole.amat[i-1]+mole.num_comole_calc;
00282 }
00283 mole.b = (double *)MALLOC((size_t)mole.num_comole_calc*sizeof(double));
00284 mole.c = (double **)MALLOC((size_t)mole.num_comole_tot*sizeof(double *));
00285 mole.c[0] = (double *)MALLOC((size_t)mole.num_comole_tot*
00286 mole.num_comole_calc*sizeof(double));
00287 for(i=1;i<mole.num_comole_tot;i++)
00288 {
00289 mole.c[i] = mole.c[i-1]+mole.num_comole_calc;
00290 }
00291 ipiv = (int32 *)MALLOC((size_t)mole.num_comole_calc*sizeof(int32));
00292
00293 tot_ion = (realnum *)MALLOC((size_t)mole.num_elements*sizeof(realnum));
00294
00295 return;
00296
00297 }
00298
00299
00300
00301 STATIC void newelement(const char label[], int ipion, int priority)
00302 {
00303 struct chem_element_s *element;
00304 data_u *p;
00305 int exists;
00306
00307 DEBUG_ENTRY("newelement()");
00308
00309 element = (struct chem_element_s *) MALLOC(sizeof(struct chem_element_s));
00310 element->ipCl = ipion;
00311 element->ipZ = priority;
00312 ASSERT ((size_t)strlen(label) < sizeof(element->chName));
00313 strncpy(element->chName,label,sizeof(element->chName));
00314 element->ipMl = element->ipMlP = -1;
00315 p = addentry(element->chName,0,mole_priv.elemtab,&exists);
00316 p->p = (void *) element;
00317 element_list[ipion] = element;
00318 }
00319
00320 STATIC struct molecule *newspecies(const char label[7], enum spectype type,
00321 enum molstate state, realnum *location, double frac0)
00322 {
00323 int exists;
00324 data_u *p;
00325 char mylab[7], thisel[CHARS_ELEMENT], *s;
00326 long int i, n, nnuc, nelem, nel;
00327 struct molecule *mol;
00328 struct chem_element_s *el, *maxel;
00329
00330 DEBUG_ENTRY("newspecies()");
00331
00332 mol = (struct molecule *) MALLOC(sizeof(struct molecule));
00333
00334 mole.num_comole_tot++;
00335 if(state == ACTIVE)
00336 mole.num_comole_calc++;
00337
00338 mol->pdr_mole_co = (realnum) frac0;
00339 for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem )
00340 {
00341 mol->nElem[nelem] = 0;
00342 }
00343 mol->co_save = 0.;
00344
00345 strncpy(mol->label,label,sizeof(mol->label));
00346 s = strchr_s(mol->label,' ');
00347 if(s)
00348 *s = '\0';
00349
00350
00351 p = addentry(mol->label,0,mole_priv.spectab,&exists);
00352 p->p = (void *) mol;
00353
00354 if(state == ACTIVE)
00355 {
00356 mol->active = 1;
00357 }
00358 else
00359 {
00360 mol->active = 0;
00361 }
00362 mol->location = location;
00363 mol->n_nuclei = 0;
00364
00365 if(type == MOLECULE)
00366 {
00367
00368 strncpy(mylab,label,7);
00369
00370
00371 mol->nElec = mol->Excit = 0;
00372
00373
00374 s = strpbrk(mylab,"*");
00375 if(s)
00376 {
00377 mol->Excit = 1;
00378 *s = '\0';
00379 }
00380
00381
00382 s = strpbrk(mylab,"+-");
00383 if(s)
00384 {
00385 if(isdigit(*(s+1)))
00386 n = atoi(s+1);
00387 else
00388 n = 1;
00389 if(*s == '+')
00390 mol->nElec = n;
00391 else
00392 mol->nElec = -n;
00393 *s = '\0';
00394 }
00395
00396 s = strstr_s(mylab,"grn");
00397 if(s)
00398 {
00399 mol->lgGas_Phase = false;
00400 *s = '\0';
00401 }
00402 else
00403 {
00404 mol->lgGas_Phase = true;
00405 }
00406
00407
00408 nnuc = 0;
00409 i = 0;
00410 maxel = NULL;
00411 while (mylab[i] != '\0' && mylab[i] != ' ' && mylab[i] != '*')
00412 {
00413 nel = 0;
00414
00415
00416 if( mylab[i] == '^' )
00417 {
00418 long A = 0;
00419 thisel[nel++] = mylab[i++];
00420 do {
00421 A = 10*A+(long int)(mylab[i]-'0');
00422 thisel[nel++] = mylab[i++];
00423 } while (isdigit(mylab[i]));
00424 }
00425
00426
00427 thisel[nel++] = mylab[i++];
00428 if(islower(mylab[i]))
00429 {
00430 thisel[nel++] = mylab[i++];
00431 }
00432 thisel[nel] = '\0';
00433
00434 el = findelement(thisel);
00435 if(el == NULL)
00436 {
00437 fprintf(stderr,"Did not recognize element at %s in \"%s \"[%ld]\n",mylab+i,label,i);
00438 cdEXIT( EXIT_FAILURE );
00439 }
00440
00441
00442
00443 if(el->ipZ != 0 && (maxel == NULL || el->ipZ > maxel->ipZ))
00444 {
00445 maxel = el;
00446 }
00447
00448 if(isdigit(mylab[i]))
00449 {
00450 n = 0;
00451 do {
00452 n = 10*n+(long int)(mylab[i]-'0');
00453 i++;
00454 } while (isdigit(mylab[i]));
00455 }
00456 else
00457 {
00458 n = 1;
00459 }
00460 mol->nElem[el->ipCl] += n;
00461 nnuc += n;
00462 }
00463
00464 mol->n_nuclei = nnuc;
00465
00466 ASSERT(mol->active == 0 || maxel != NULL );
00467
00468 if(maxel != NULL)
00469 mol->nelem_hevmol = maxel->ipCl;
00470 else
00471 mol->nelem_hevmol = -1;
00472
00473 mol->mole_mass = 0.;
00474 for(i=0;i<LIMELM;++i)
00475 {
00476 mol->mole_mass +=
00477 mol->nElem[i]*dense.AtomicWeight[i]*(realnum)ATOMIC_MASS_UNIT;
00478 }
00479
00480
00481 if(mol->n_nuclei == 1 && mol->active && mol->nElec == 0)
00482 {
00483 mole.num_elements++;
00484 mole.lgElem_in_chemistry[mol->nelem_hevmol] = true;
00485 }
00486 }
00487 return mol;
00488 }
00489 STATIC int isactive(data_u *dat)
00490 {
00491
00492 DEBUG_ENTRY("isactive()");
00493
00494 struct molecule *p = (struct molecule *) (dat->p);
00495
00496 return p->active;
00497 }
00498 STATIC int ispassive(data_u *dat)
00499 {
00500 return !((struct molecule *) dat->p)->active;
00501 }
00502 STATIC int isCOnet(data_u *dat)
00503 {
00504 return (((struct chem_element_s *) dat->p)->ipMl != -1);
00505 }
00506
00507 struct molecule *findspecies(const char buf[])
00508 {
00509 DEBUG_ENTRY("findspecies()");
00510
00511
00512 string s;
00513 for (const char *pb = buf; *pb && *pb != ' '; ++pb)
00514 {
00515 s += *pb;
00516 }
00517
00518 data_u *p = lookup(s.c_str(),0,mole_priv.spectab);
00519
00520 if(p)
00521 return (struct molecule *) p->p;
00522 else
00523 return &null_mole;
00524 }
00525
00526 STATIC struct chem_element_s *findelement(const char buf[])
00527 {
00528 data_u *p;
00529
00530 DEBUG_ENTRY("findelement()");
00531
00532 p = lookup(buf,0,mole_priv.elemtab);
00533
00534 if(p)
00535 return (struct chem_element_s *) p->p;
00536 else
00537 return NULL;
00538 }
00539
00540 struct COmole_rate_s *CO_findrate_s(const char buf[])
00541 {
00542 data_u *p;
00543
00544 DEBUG_ENTRY("CO_findrate_s()");
00545
00546 p = lookup(buf,0,mole_priv.reactab);
00547
00548 if(p)
00549 return (struct COmole_rate_s *) p->p;
00550 else
00551 return NULL;
00552 }
00553 double CO_findrk(const char buf[])
00554 {
00555 struct COmole_rate_s *rate;
00556
00557 DEBUG_ENTRY("CO_findrk()");
00558
00559 rate = CO_findrate_s(buf);
00560
00561 if(!rate)
00562 return 0.;
00563
00564
00565 ASSERT( !isnan( rate->rk ) );
00566
00567 return rate->rk;
00568 }
00569 double CO_findrate(const char buf[])
00570 {
00571 struct COmole_rate_s *rate;
00572 double ret;
00573 int i;
00574
00575 DEBUG_ENTRY("CO_findrate()");
00576
00577 rate = CO_findrate_s(buf);
00578 if(!rate)
00579 {
00580 return 0.;
00581 }
00582
00583 ret = rate->rk;
00584 for(i=0;i<rate->nrates;i++)
00585 ret *= rate->rate_species[i]->hevmol;
00586 return ret;
00587 }
00588
00589 void CO_update_species_cache(void)
00590 {
00591 int i;
00592
00593 DEBUG_ENTRY("CO_update_species_cache()");
00594
00595 dense.eden_f = (realnum)dense.eden;
00596
00597 for(i=0;i<mole.num_comole_tot;i++)
00598 {
00599 if(COmole[i]->location)
00600 {
00601 COmole[i]->hevmol = *(COmole[i]->location);
00602
00603
00604 ASSERT( !isnan( COmole[i]->hevmol ) );
00605
00606 ASSERT( COmole[i]->hevmol < MAX_DENSITY );
00607 }
00608 }
00609 }
00610
00611
00612
00613
00614 double CO_sink_rate(const char chSpecies[7])
00615 {
00616 int ipthis, i, n, nt;
00617 double ratev, ratevi;
00618 struct COmole_rate_s *rate;
00619 struct molecule *sp;
00620
00621 DEBUG_ENTRY("CO_sink_rate()");
00622
00623 sp = findspecies(chSpecies);
00624 ratev = 0;
00625 nt = coreactions.n;
00626
00627 for(n=0;n<nt;n++)
00628 {
00629 rate = coreactions.list[n];
00630 ipthis = -1;
00631 for(i=0;i<rate->nrates && ipthis == -1;i++)
00632 {
00633 if(rate->rate_species[i] == sp)
00634 {
00635 ipthis = i;
00636 }
00637 }
00638 if(ipthis != -1) {
00639 ratevi = rate->rk;
00640 for(i=0;i<rate->nrates;i++)
00641 {
00642 if(i!=ipthis)
00643 {
00644 ratevi *= rate->rate_species[i]->hevmol;
00645 }
00646 }
00647 ratev += ratevi;
00648 }
00649 }
00650
00651 return ratev;
00652 }
00653 #ifdef PRINTREACT
00654 STATIC void printreact(struct COmole_rate_s *rate)
00655 {
00656 int i;
00657
00658 DEBUG_ENTRY("printreact()");
00659
00660 for(i=0;i<rate->nreactants;i++) {
00661 fprintf(stderr,"%s,",rate->reactants[i]->label);
00662 }
00663 fprintf(stderr,"=>");
00664 for(i=0;i<rate->nproducts;i++) {
00665 fprintf(stderr,"%s,",rate->products[i]->label);
00666 }
00667 fprintf(stderr,"\n");
00668
00669 }
00670 #endif
00671 void CO_update_rks(void)
00672 {
00673 int n, nt;
00674 struct COmole_rate_s *rate;
00675
00676 DEBUG_ENTRY("CO_update_rks()");
00677
00678 nt = coreactions.n;
00679
00680 for(n=0;n<nt;n++)
00681 {
00682 rate = coreactions.list[n];
00683 if(rate->fun != NULL) {
00684 rate->rk = rate->a*rate->fun(rate);
00685 }
00686 }
00687 }
00688
00689 double CO_dissoc_rate(const char chSpecies[7])
00690 {
00691 int ipthis, i, n, nt;
00692 double ratev, ratevi;
00693 struct COmole_rate_s *rate;
00694 struct molecule *sp;
00695
00696 DEBUG_ENTRY("CO_dissoc_rate()");
00697
00698 sp = findspecies(chSpecies);
00699 ratev = 0;
00700 nt = coreactions.n;
00701
00702 for(n=0;n<nt;n++)
00703 {
00704 rate = coreactions.list[n];
00705 if(rate->photon == -1)
00706 {
00707 ipthis = 0;
00708 for(i=0;i<rate->nproducts;i++)
00709 {
00710 if(rate->products[i] == sp)
00711 {
00712 ipthis++;
00713 }
00714 }
00715 if(ipthis)
00716 {
00717 ratevi = rate->rk;
00718 for(i=0;i<rate->nrates;i++)
00719 {
00720 ratevi *= rate->rate_species[i]->hevmol;
00721 }
00722 ratev += ipthis*ratevi;
00723 }
00724 }
00725 }
00726
00727 return ratev;
00728 }
00729 double CO_source_rate(const char chSpecies[7])
00730 {
00731 int ipthis, i, n, nt;
00732 double ratev, ratevi;
00733 struct COmole_rate_s *rate;
00734 struct molecule *sp;
00735
00736 DEBUG_ENTRY("CO_source_rate()");
00737
00738 sp = findspecies(chSpecies);
00739 ratev = 0;
00740 nt = coreactions.n;
00741
00742 for(n=0;n<nt;n++)
00743 {
00744 rate = coreactions.list[n];
00745 ipthis = 0;
00746 for(i=0;i<rate->nproducts;i++)
00747 {
00748 if(rate->products[i] == sp)
00749 {
00750 ipthis++;
00751 }
00752 }
00753 if(ipthis) {
00754 ratevi = rate->rk;
00755 for(i=0;i<rate->nrates;i++)
00756 {
00757 ratevi *= rate->rate_species[i]->hevmol;
00758 }
00759 ratev += ipthis*ratevi;
00760 }
00761 }
00762
00763 return ratev;
00764 }
00765
00766
00767 void CO_punch_mol(FILE *punit, const char chSpecies[], char header[], double depth)
00768 {
00769 int n, nt, i, ipthis;
00770 double ratevi;
00771 struct COmole_rate_s *rate;
00772 struct molecule *sp;
00773 char *s;
00774
00775 DEBUG_ENTRY("CO_punch_mol()");
00776
00777 s = header;
00778
00779 sp = findspecies(chSpecies);
00780 if(punit == NULL)
00781 {
00782 if (sp == &null_mole)
00783 {
00784 fprintf(ioQQQ,"Could not find species '%s' to save chemistry rates.\nSorry.\n",
00785 chSpecies);
00786 cdEXIT(EXIT_FAILURE);
00787 }
00788 sprintf (s,"#Depth");
00789 s += strlen(s);
00790 }
00791 else
00792 {
00793 fprintf (punit,"%.5e",depth);
00794 }
00795
00796 nt = coreactions.n;
00797 for(n=0;n<nt;n++) {
00798 rate = coreactions.list[n];
00799 ipthis = 0;
00800 for(i=0;i<rate->nrates;i++)
00801 {
00802 if(rate->rate_species[i] == sp)
00803 {
00804 ipthis++;
00805 }
00806 }
00807 for(i=0;i<rate->nproducts;i++)
00808 {
00809 if(rate->products[i] == sp)
00810 {
00811 ipthis++;
00812 }
00813 }
00814 if(ipthis)
00815 {
00816 if(punit == NULL)
00817 {
00818 sprintf(s,"\t%s",rate->label);
00819 s += strlen(s);
00820 }
00821 else
00822 {
00823 ratevi = rate->rk;
00824 for(i=0;i<rate->nrates;i++)
00825 {
00826 ratevi *= rate->rate_species[i]->hevmol;
00827 }
00828 fprintf(punit,"\t%.3e",ratevi);
00829 }
00830 }
00831 }
00832 if(punit == NULL)
00833 {
00834 sprintf(s,"\n");
00835 }
00836 else
00837 {
00838 fprintf(punit,"\n");
00839 }
00840 }
00841
00842 void CO_zero(void)
00843 {
00844 long int i;
00845
00846 static bool lgFirstCall = true;
00847 static long int num_comole_calc_MALLOC=-1;
00848
00849 DEBUG_ENTRY("CO_zero()");
00850
00851 if( lgFirstCall )
00852 {
00853
00854 timesc.AgeCOMoleDest = (double*)MALLOC( mole.num_comole_calc*sizeof(double) );
00855
00856 lgFirstCall = false;
00857 num_comole_calc_MALLOC = mole.num_comole_calc;
00858 }
00859 else if( mole.num_comole_calc>num_comole_calc_MALLOC )
00860 {
00861
00862
00863 fprintf(ioQQQ,"DISASTER - the number of species in the CO network has increased. This is not allowed.\n");
00864 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");
00865 fprintf(ioQQQ,"Sorry.\n");
00866 cdEXIT(EXIT_FAILURE);
00867 }
00868
00869 for( i=0; i < mole.num_comole_calc; i++ )
00870 {
00871 timesc.AgeCOMoleDest[i] = 0.;
00872 }
00873
00874
00875 for( i=0; i<mole.num_comole_calc; ++i )
00876 COmole[i]->xMoleFracMax = 0.;
00877
00878
00879 for( i=0; i < mole.num_comole_tot; i++ )
00880 {
00881 COmole[i]->hevmol = 0.;
00882 COmole[i]->hevcol = 0.;
00883 }
00884 }