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
00021 enum spectype {MOLECULE, OTHER};
00022 enum molstate {ACTIVE, PASSIVE};
00023
00024 static void newelement(const char label[], int ipion,
00025 int priority);
00026 static struct molecule *newspecies(const char label[7], enum spectype type,
00027 enum molstate state, float *location, double frac0);
00028 static struct chem_element_s *findelement(const char buf[]);
00029 static int isactive(data_u *dat);
00030 static int ispassive(data_u *dat);
00031 static int isCOnet(data_u *dat);
00032
00033 struct chem_element_s **chem_element;
00034
00035
00036 struct chem_element_s *element_list[LIMELM];
00037 int32 *ipiv;
00038 float *tot_ion;
00039
00040 struct mole_priv_s mole_priv;
00041 struct molecule null_mole;
00042
00043
00044
00045 void CO_Init(void)
00046 {
00047 long int i,
00048 nelem;
00049 struct molecule *sp;
00050 static bool lgCO_Init_called=false;
00051
00052 DEBUG_ENTRY( "CO_Init()" );
00053
00054
00055
00056 co.co_nzone = -2;
00057 co.iteration_co = -2;
00058
00059
00060
00061
00062 if( lgCO_Init_called )
00063 {
00064 DEBUG_EXIT( "CO_Init()" );
00065 return;
00066 }
00067
00068
00069 lgCO_Init_called = true;
00070
00071 mole_priv.spectab = newhash(NULL);
00072 mole_priv.reactab = newhash(NULL);
00073 mole_priv.elemtab = newhash(NULL);
00074
00075 for(nelem=ipHYDROGEN; nelem<LIMELM; ++nelem )
00076 {
00077 element_list[nelem] = NULL;
00078 mole.lgElem_in_chemistry[nelem] = false;
00079 }
00080
00081
00082
00083
00084
00085
00086
00087 newelement("C" ,ipCARBON,3);
00088 newelement("O" ,ipOXYGEN,2);
00089 newelement("Si",ipSILICON,5);
00090 newelement("N" ,ipNITROGEN,4);
00091 newelement("S" ,ipSULPHUR,6);
00092 newelement("Cl",ipCHLORINE,7);
00093 newelement("H" ,ipHYDROGEN,1);
00094 newelement("He",ipHELIUM,0);
00095 newelement("Mg",ipMAGNESIUM,0);
00096 newelement("Fe",ipIRON,0);
00097
00098
00099
00100
00101
00102
00103
00104
00105
00106
00107 mole.num_comole_calc = mole.num_comole_tot = mole.num_elements = 0;
00108 null_mole.hevmol = null_mole.hevcol = 0.;
00109 null_mole.index = -1;
00110
00111 sp = newspecies("C ",MOLECULE,ACTIVE,NULL,1.00e+00);
00112 sp = newspecies("O ",MOLECULE,ACTIVE,NULL,1.13e+05);
00113 sp = newspecies("Si ",MOLECULE,ACTIVE,NULL,3.56e-04);
00114 sp = newspecies("N ",MOLECULE,ACTIVE,NULL,2.25e+00);
00115 sp = newspecies("S ",MOLECULE,ACTIVE,NULL,6.90e-03);
00116 sp = newspecies("Cl ",MOLECULE,ACTIVE,NULL,6.90e-03);
00117 sp = newspecies("C+ ",MOLECULE,ACTIVE,NULL,6.79e+04);
00118 sp = newspecies("O+ ",MOLECULE,ACTIVE,NULL,1.58e+01);
00119 sp = newspecies("Si+ ",MOLECULE,ACTIVE,NULL,1.79e+02);
00120 sp = newspecies("N+ ",MOLECULE,ACTIVE,NULL,2.62e-10);
00121 sp = newspecies("S+ ",MOLECULE,ACTIVE,NULL,1.79e+03);
00122 sp = newspecies("Cl+ ",MOLECULE,ACTIVE,NULL,3.71e-09);
00123 sp = newspecies("CH ",MOLECULE,ACTIVE,NULL,1.10e-06);
00124 sp = newspecies("CH+ ",MOLECULE,ACTIVE,NULL,6.99e-05);
00125 sp = newspecies("OH ",MOLECULE,ACTIVE,NULL,2.15e-03);
00126 sp = newspecies("OH+ ",MOLECULE,ACTIVE,NULL,2.45e-04);
00127 sp = newspecies("O2 ",MOLECULE,ACTIVE,NULL,1.91e-07);
00128 sp = newspecies("CO ",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 = (float)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 = (float *)MALLOC((size_t)mole.num_elements*sizeof(float));
00294
00295 DEBUG_EXIT( "CO_Init()" );
00296
00297 return;
00298
00299 }
00300
00301
00302
00303
00304 static void newelement(const char label[], int ipion, int priority)
00305 {
00306 struct chem_element_s *element;
00307 data_u *p;
00308 int exists;
00309
00310 DEBUG_ENTRY("newelement()");
00311
00312 element = (struct chem_element_s *) MALLOC(sizeof(struct chem_element_s));
00313 element->ipCl = ipion;
00314 element->ipZ = priority;
00315 ASSERT ((size_t)strlen(label) < sizeof(element->chName));
00316 strncpy(element->chName,label,sizeof(element->chName));
00317 element->ipMl = element->ipMlP = -1;
00318 p = addentry(element->chName,0,mole_priv.elemtab,&exists);
00319 p->p = (void *) element;
00320 element_list[ipion] = element;
00321
00322 DEBUG_EXIT("newelement()");
00323 }
00324
00325 #include <string.h>
00326 #include <ctype.h>
00327 static struct molecule *newspecies(const char label[7], enum spectype type,
00328 enum molstate state, float *location, double frac0)
00329 {
00330 int exists;
00331 data_u *p;
00332 char mylab[7], thisel[3], *s;
00333 long int i, n, nnuc, nelem, nel;
00334 struct molecule *mol;
00335 struct chem_element_s *el, *maxel;
00336
00337 DEBUG_ENTRY("newspecies()");
00338
00339 mol = (struct molecule *) MALLOC(sizeof(struct molecule));
00340
00341 mole.num_comole_tot++;
00342 if(state == ACTIVE)
00343 mole.num_comole_calc++;
00344
00345 mol->pdr_mole_co = (float) frac0;
00346 for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem )
00347 {
00348 mol->nElem[nelem] = 0;
00349 }
00350 mol->co_save = 0.;
00351
00352 strncpy(mol->label,label,sizeof(mol->label));
00353 s = strchr(mol->label,' ');
00354 if(s)
00355 *s = '\0';
00356
00357
00358 p = addentry(mol->label,0,mole_priv.spectab,&exists);
00359 p->p = (void *) mol;
00360
00361 if(state == ACTIVE)
00362 {
00363 mol->active = 1;
00364 }
00365 else
00366 {
00367 mol->active = 0;
00368 }
00369 mol->location = location;
00370 mol->n_nuclei = 0;
00371
00372 if(type == MOLECULE)
00373 {
00374
00375 strncpy(mylab,label,7);
00376
00377
00378 mol->nElec = mol->Excit = 0;
00379
00380
00381 s = strpbrk(mylab,"*");
00382 if(s)
00383 {
00384 mol->Excit = 1;
00385 *s = '\0';
00386 }
00387
00388
00389 s = strpbrk(mylab,"+-");
00390 if(s)
00391 {
00392 if(isdigit(*(s+1)))
00393 n = atoi(s+1);
00394 else
00395 n = 1;
00396 if(*s == '+')
00397 mol->nElec = n;
00398 else
00399 mol->nElec = -n;
00400 *s = '\0';
00401 }
00402
00403 s = strstr(mylab,"grn");
00404 if(s)
00405 {
00406 mol->lgGas_Phase = false;
00407 *s = '\0';
00408 }
00409 else
00410 {
00411 mol->lgGas_Phase = true;
00412 }
00413
00414
00415 nnuc = 0;
00416 i = 0;
00417 maxel = NULL;
00418 while (mylab[i] != '\0' && mylab[i] != ' ' && mylab[i] != '*')
00419 {
00420
00421 nel = 0;
00422 thisel[nel++] = mylab[i++];
00423 if(islower(mylab[i]))
00424 {
00425 thisel[nel++] = mylab[i++];
00426 }
00427 thisel[nel] = '\0';
00428
00429 el = findelement(thisel);
00430 if(el == NULL)
00431 {
00432 fprintf(stderr,"Did not recognize element at %s in \"%s \"[%ld]\n",mylab+i,label,i);
00433 exit(-1);
00434 }
00435
00436
00437
00438 if(el->ipZ != 0 && (maxel == NULL || el->ipZ > maxel->ipZ))
00439 {
00440 maxel = el;
00441 }
00442
00443 if(isdigit(mylab[i]))
00444 {
00445 n = 0;
00446 do {
00447 n = 10*n+(long int)(mylab[i]-'0');
00448 i++;
00449 } while (isdigit(mylab[i]));
00450 }
00451 else
00452 {
00453 n = 1;
00454 }
00455 mol->nElem[el->ipCl] += n;
00456 nnuc += n;
00457 }
00458
00459 mol->n_nuclei = nnuc;
00460
00461 ASSERT(mol->active == 0 || maxel != NULL );
00462
00463 if(maxel != NULL)
00464 mol->nelem_hevmol = maxel->ipCl;
00465 else
00466 mol->nelem_hevmol = -1;
00467
00468 mol->mole_mass = 0.;
00469 for(i=0;i<LIMELM;++i)
00470 {
00471 mol->mole_mass +=
00472 mol->nElem[i]*dense.AtomicWeight[i]*(float)ATOMIC_MASS_UNIT;
00473 }
00474
00475
00476 if(mol->n_nuclei == 1 && mol->active && mol->nElec == 0)
00477 {
00478 mole.num_elements++;
00479 mole.lgElem_in_chemistry[mol->nelem_hevmol] = true;
00480 }
00481 }
00482
00483 DEBUG_EXIT("newspecies()");
00484 return mol;
00485 }
00486 static int isactive(data_u *dat)
00487 {
00488
00489 DEBUG_ENTRY("isactive()");
00490
00491 struct molecule *p = (struct molecule *) (dat->p);
00492
00493 DEBUG_EXIT("isactive()");
00494
00495 return p->active;
00496 }
00497 static int ispassive(data_u *dat)
00498 {
00499 return !((struct molecule *) dat->p)->active;
00500 }
00501 static int isCOnet(data_u *dat)
00502 {
00503 return (((struct chem_element_s *) dat->p)->ipMl != -1);
00504 }
00505
00506 struct molecule *findspecies(const char buf[])
00507 {
00508 data_u *p;
00509
00510 DEBUG_ENTRY("findspecies()");
00511
00512 p = lookup(buf,0,mole_priv.spectab);
00513
00514 DEBUG_EXIT("findspecies()");
00515
00516 if(p)
00517 return (struct molecule *) p->p;
00518 else
00519 return &null_mole;
00520 }
00521
00522 static struct chem_element_s *findelement(const char buf[])
00523 {
00524 data_u *p;
00525
00526 DEBUG_ENTRY("findelement()");
00527
00528 p = lookup(buf,0,mole_priv.elemtab);
00529
00530 DEBUG_EXIT("findelement()");
00531
00532 if(p)
00533 return (struct chem_element_s *) p->p;
00534 else
00535 return NULL;
00536 }
00537
00538 struct COmole_rate_s *CO_findrate_s(const char buf[])
00539 {
00540 data_u *p;
00541
00542 DEBUG_ENTRY("CO_findrate_s()");
00543
00544 p = lookup(buf,0,mole_priv.reactab);
00545
00546 DEBUG_EXIT("CO_findrate_s()");
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 DEBUG_EXIT("CO_findrk()");
00562
00563 if(!rate)
00564 return 0.;
00565
00566
00567 ASSERT( !isnan( rate->rk ) );
00568
00569 return rate->rk;
00570 }
00571 double CO_findrate(const char buf[])
00572 {
00573 struct COmole_rate_s *rate;
00574 double ret;
00575 int i;
00576
00577 DEBUG_ENTRY("CO_findrate()");
00578
00579 rate = CO_findrate_s(buf);
00580 if(!rate)
00581 {
00582 DEBUG_EXIT("CO_findrate()");
00583 return 0.;
00584 }
00585
00586 ret = rate->rk;
00587 for(i=0;i<rate->nrates;i++)
00588 ret *= rate->rate_species[i]->hevmol;
00589
00590 DEBUG_EXIT("CO_findrate()");
00591 return ret;
00592 }
00593
00594 void CO_update_species_cache(void)
00595 {
00596 int i;
00597
00598 DEBUG_ENTRY("CO_update_species_cache()");
00599
00600 dense.eden_f = (float)dense.eden;
00601
00602 for(i=0;i<mole.num_comole_tot;i++)
00603 {
00604 if(COmole[i]->location)
00605 {
00606 COmole[i]->hevmol = *(COmole[i]->location);
00607
00608
00609 ASSERT( !isnan( COmole[i]->hevmol ) );
00610 }
00611 }
00612
00613 DEBUG_EXIT("CO_update_species_cache()");
00614 }
00615
00616
00617
00618
00619 double CO_sink_rate(const char species[7])
00620 {
00621 int ipthis, i, n, nt;
00622 double ratev, ratevi;
00623 struct COmole_rate_s *rate;
00624 struct molecule *sp;
00625
00626 DEBUG_ENTRY("CO_sink_rate()");
00627
00628 sp = findspecies(species);
00629 ratev = 0;
00630 nt = coreactions.n;
00631
00632 for(n=0;n<nt;n++)
00633 {
00634 rate = coreactions.list[n];
00635 ipthis = -1;
00636 for(i=0;i<rate->nrates && ipthis == -1;i++)
00637 {
00638 if(rate->rate_species[i] == sp)
00639 {
00640 ipthis = i;
00641 }
00642 }
00643 if(ipthis != -1) {
00644 ratevi = rate->rk;
00645 for(i=0;i<rate->nrates;i++)
00646 {
00647 if(i!=ipthis)
00648 {
00649 ratevi *= rate->rate_species[i]->hevmol;
00650 }
00651 }
00652 ratev += ratevi;
00653 }
00654 }
00655
00656 DEBUG_EXIT("CO_sink_rate()");
00657
00658 return ratev;
00659 }
00660 #ifdef PRINTREACT
00661 static void printreact(struct COmole_rate_s *rate)
00662 {
00663 int i;
00664
00665 DEBUG_ENTRY("printreact()");
00666
00667 for(i=0;i<rate->nreactants;i++) {
00668 fprintf(stderr,"%s,",rate->reactants[i]->label);
00669 }
00670 fprintf(stderr,"=>");
00671 for(i=0;i<rate->nproducts;i++) {
00672 fprintf(stderr,"%s,",rate->products[i]->label);
00673 }
00674 fprintf(stderr,"\n");
00675
00676 DEBUG_EXIT("printreact()");
00677
00678 }
00679 #endif
00680 void CO_update_rks(void)
00681 {
00682 int n, nt;
00683 struct COmole_rate_s *rate;
00684
00685 DEBUG_ENTRY("CO_update_rks()");
00686
00687 nt = coreactions.n;
00688
00689 for(n=0;n<nt;n++)
00690 {
00691 rate = coreactions.list[n];
00692 if(rate->fun != NULL) {
00693 rate->rk = rate->a*rate->fun(rate);
00694 }
00695 }
00696
00697 DEBUG_EXIT("CO_update_rks()");
00698 }
00699
00700 double CO_dissoc_rate(const char species[7])
00701 {
00702 int ipthis, i, n, nt;
00703 double ratev, ratevi;
00704 struct COmole_rate_s *rate;
00705 struct molecule *sp;
00706
00707 DEBUG_ENTRY("CO_dissoc_rate()");
00708
00709 sp = findspecies(species);
00710 ratev = 0;
00711 nt = coreactions.n;
00712
00713 for(n=0;n<nt;n++)
00714 {
00715 rate = coreactions.list[n];
00716 if(rate->photon == -1)
00717 {
00718 ipthis = 0;
00719 for(i=0;i<rate->nproducts;i++)
00720 {
00721 if(rate->products[i] == sp)
00722 {
00723 ipthis++;
00724 }
00725 }
00726 if(ipthis)
00727 {
00728 ratevi = rate->rk;
00729 for(i=0;i<rate->nrates;i++)
00730 {
00731 ratevi *= rate->rate_species[i]->hevmol;
00732 }
00733 ratev += ipthis*ratevi;
00734 }
00735 }
00736 }
00737
00738 DEBUG_EXIT("CO_dissoc_rate()");
00739
00740 return ratev;
00741 }
00742 double CO_source_rate(const char species[7])
00743 {
00744 int ipthis, i, n, nt;
00745 double ratev, ratevi;
00746 struct COmole_rate_s *rate;
00747 struct molecule *sp;
00748
00749 DEBUG_ENTRY("CO_source_rate()");
00750
00751 sp = findspecies(species);
00752 ratev = 0;
00753 nt = coreactions.n;
00754
00755 for(n=0;n<nt;n++)
00756 {
00757 rate = coreactions.list[n];
00758 ipthis = 0;
00759 for(i=0;i<rate->nproducts;i++)
00760 {
00761 if(rate->products[i] == sp)
00762 {
00763 ipthis++;
00764 }
00765 }
00766 if(ipthis) {
00767 ratevi = rate->rk;
00768 for(i=0;i<rate->nrates;i++)
00769 {
00770 ratevi *= rate->rate_species[i]->hevmol;
00771 }
00772 ratev += ipthis*ratevi;
00773 }
00774 }
00775
00776 DEBUG_EXIT("CO_source_rate()");
00777
00778 return ratev;
00779 }
00780
00781
00782 void CO_punch_mol(FILE *punit, const char species[], char header[], double depth)
00783 {
00784 int n, nt, i, ipthis;
00785 double ratevi;
00786 struct COmole_rate_s *rate;
00787 struct molecule *sp;
00788 char *s;
00789
00790 DEBUG_ENTRY("CO_punch_mol()");
00791
00792 s = header;
00793
00794 sp = findspecies(species);
00795 if(punit == NULL)
00796 {
00797 sprintf (s,"#Depth");
00798 s += strlen(s);
00799 }
00800 else
00801 {
00802 fprintf (punit,"%.5e",depth);
00803 }
00804
00805 nt = coreactions.n;
00806 for(n=0;n<nt;n++) {
00807 rate = coreactions.list[n];
00808 ipthis = 0;
00809 for(i=0;i<rate->nrates;i++)
00810 {
00811 if(rate->rate_species[i] == sp)
00812 {
00813 ipthis++;
00814 }
00815 }
00816 for(i=0;i<rate->nproducts;i++)
00817 {
00818 if(rate->products[i] == sp)
00819 {
00820 ipthis++;
00821 }
00822 }
00823 if(ipthis)
00824 {
00825 if(punit == NULL)
00826 {
00827 sprintf(s,"\t%s",rate->label);
00828 s += strlen(s);
00829 }
00830 else
00831 {
00832 ratevi = rate->rk;
00833 for(i=0;i<rate->nrates;i++)
00834 {
00835 ratevi *= rate->rate_species[i]->hevmol;
00836 }
00837 fprintf(punit,"\t%.3e",ratevi);
00838 }
00839 }
00840 }
00841 if(punit == NULL)
00842 {
00843 sprintf(s,"\n");
00844 }
00845 else
00846 {
00847 fprintf(punit,"\n");
00848 }
00849
00850 DEBUG_EXIT("CO_punch_mol()");
00851 }
00852
00853 void CO_zero(void)
00854 {
00855 long int i;
00856
00857 static bool lgFirstCall = true;
00858 static long int num_comole_calc_MALLOC=-1;
00859
00860 DEBUG_ENTRY("CO_zero()");
00861
00862 if( lgFirstCall )
00863 {
00864
00865 timesc.AgeCOMoleDest = (double*)MALLOC( mole.num_comole_calc*sizeof(double) );
00866
00867 lgFirstCall = false;
00868 num_comole_calc_MALLOC = mole.num_comole_calc;
00869 }
00870 else if( mole.num_comole_calc>num_comole_calc_MALLOC )
00871 {
00872
00873
00874 fprintf(ioQQQ,"DISASTER - the number of species in the CO network has increased. This is not allowed.\n");
00875 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");
00876 fprintf(ioQQQ,"Sorry.\n");
00877 cdEXIT(EXIT_FAILURE);;
00878 }
00879
00880 for( i=0; i < mole.num_comole_calc; i++ )
00881 {
00882 timesc.AgeCOMoleDest[i] = 0.;
00883 }
00884
00885
00886 for( i=0; i<mole.num_comole_calc; ++i )
00887 COmole[i]->xMoleFracMax = 0.;
00888
00889
00890 for( i=0; i < mole.num_comole_tot; i++ )
00891 {
00892 COmole[i]->hevmol = 0.;
00893 COmole[i]->hevcol = 0.;
00894 }
00895
00896 DEBUG_EXIT("CO_zero()");
00897 }
00898