cloudy trunk
|
00001 /* This file is part of Cloudy and is copyright (C)1978-2008 by Gary J. Ferland and 00002 * others. For conditions of distribution and use see copyright notice in license.txt */ 00003 /*CO_Init called from cdInit to initialize co routines */ 00004 /*CO_update_chem_rates update rate coefficients, only temp part - in mole_co_etc.c */ 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 /*lint -e778 constant expression evaluates to 0 in operation '-' */ 00016 00017 /*CO_update_chem_rates update rate coefficients, only temp part - in mole_co_etc.c 00018 * called in conv_base before any chemistry or ionization is done */ 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 /* List of element structures indexed by atom index 00034 -- could use (element_list[nelem] != NULL) for mole.lgElem_in_chemistry[nelem] */ 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 /*CO_Init called from cdInit to initialize CO routines */ 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 /* these tell the molecular solver what zone and iteration it has 00054 * been evaluated on */ 00055 co.co_nzone = -2; 00056 co.iteration_co = -2; 00057 00058 /* prevent memory leaks */ 00059 /* \todo this is a temporary fix for PR14. We should improve the overall design 00060 * of this code to prevent valid pointers being overwritten in a second call to CO_Init */ 00061 if( lgCO_Init_called ) 00062 { 00063 return; 00064 } 00065 00066 /* say that we have been called */ 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 /* set up concordance of elemental species to external Cloudy indices */ 00080 /* final number is 'element priority':- 00081 rjrw 2006 Aug 11: Nick Abel explains this order is roughly gas phase 00082 abundance -- the encoding of atomic species is from the last 00083 elements of "enum molecule_codes" in mole.h. Should this make any 00084 difference? In practice it does. */ 00085 newelement("C" ,ipCARBON,3); 00086 newelement("O" ,ipOXYGEN,2); 00087 newelement("Si",ipSILICON,5); 00088 newelement("N" ,ipNITROGEN,4); 00089 newelement("S" ,ipSULPHUR,6); 00090 newelement("Cl",ipCHLORINE,7); 00091 newelement("H" ,ipHYDROGEN,1); 00092 newelement("He",ipHELIUM,0); 00093 newelement("Mg",ipMAGNESIUM,0); 00094 newelement("Fe",ipIRON,0); 00095 00096 /* set up properties of molecular species -- chemical formulae, 00097 array indices, elementary components (parsed from formula), 00098 status within CO network, location of stored value external 00099 to CO network (must be floating point). */ 00100 00101 /* final arguments are relative to neutral carbon and are the solution to the first zone of the 00102 * TH85 pdr */ 00103 00104 /* Sizes of different parts of network are calculated in newspecies */ 00105 mole.num_comole_calc = mole.num_comole_tot = mole.num_elements = 0; 00106 null_mole.hevmol = null_mole.hevcol = 0.; /* Non-molecule to allow valid pointer return inactive species are accessed */ 00107 null_mole.index = -1; 00108 00109 sp = newspecies("C ",MOLECULE,ACTIVE,NULL,1.00e+00); 00110 sp = newspecies("O ",MOLECULE,ACTIVE,NULL,1.13e+05); 00111 sp = newspecies("Si ",MOLECULE,ACTIVE,NULL,3.56e-04); 00112 sp = newspecies("N ",MOLECULE,ACTIVE,NULL,2.25e+00); 00113 sp = newspecies("S ",MOLECULE,ACTIVE,NULL,6.90e-03); 00114 sp = newspecies("Cl ",MOLECULE,ACTIVE,NULL,6.90e-03); 00115 sp = newspecies("C+ ",MOLECULE,ACTIVE,NULL,6.79e+04); 00116 sp = newspecies("O+ ",MOLECULE,ACTIVE,NULL,1.58e+01); 00117 sp = newspecies("Si+ ",MOLECULE,ACTIVE,NULL,1.79e+02); 00118 sp = newspecies("N+ ",MOLECULE,ACTIVE,NULL,2.62e-10); 00119 sp = newspecies("S+ ",MOLECULE,ACTIVE,NULL,1.79e+03); 00120 sp = newspecies("Cl+ ",MOLECULE,ACTIVE,NULL,3.71e-09); 00121 sp = newspecies("CH ",MOLECULE,ACTIVE,NULL,1.10e-06); 00122 sp = newspecies("CH+ ",MOLECULE,ACTIVE,NULL,6.99e-05); 00123 sp = newspecies("OH ",MOLECULE,ACTIVE,NULL,2.15e-03); 00124 sp = newspecies("OH+ ",MOLECULE,ACTIVE,NULL,2.45e-04); 00125 sp = newspecies("O2 ",MOLECULE,ACTIVE,NULL,1.91e-07); 00126 sp = newspecies("CO ",MOLECULE,ACTIVE,NULL,4.05e-05); 00127 sp = newspecies("CO+ ",MOLECULE,ACTIVE,NULL,2.16e-06); 00128 sp = newspecies("H2O ",MOLECULE,ACTIVE,NULL,1.60e-11); 00129 sp = newspecies("H2O+ ",MOLECULE,ACTIVE,NULL,1.27e-10); 00130 sp = newspecies("O2+ ",MOLECULE,ACTIVE,NULL,1.17e-06); 00131 sp = newspecies("H3O+ ",MOLECULE,ACTIVE,NULL,3.44e-17); 00132 sp = newspecies("CH2+ ",MOLECULE,ACTIVE,NULL,3.71e-09); 00133 sp = newspecies("CH2 ",MOLECULE,ACTIVE,NULL,8.59e-13); 00134 sp = newspecies("HCO+ ",MOLECULE,ACTIVE,NULL,4.01e-10); 00135 sp = newspecies("CH3+ ",MOLECULE,ACTIVE,NULL,7.02e-16); 00136 sp = newspecies("CH3 ",MOLECULE,ACTIVE,NULL,4.59e-19); 00137 sp = newspecies("CH4 ",MOLECULE,ACTIVE,NULL,9.12e-28); 00138 sp = newspecies("CH4+ ",MOLECULE,ACTIVE,NULL,1.03e-28); 00139 sp = newspecies("CH5+ ",MOLECULE,ACTIVE,NULL,8.16e-28); 00140 sp = newspecies("SiH2+ ",MOLECULE,ACTIVE,NULL,2.07e-13); 00141 sp = newspecies("SiH ",MOLECULE,ACTIVE,NULL,1.80e-14); 00142 sp = newspecies("SiOH+ ",MOLECULE,ACTIVE,NULL,5.39e-15); 00143 sp = newspecies("SiO ",MOLECULE,ACTIVE,NULL,3.89e-14); 00144 sp = newspecies("SiO+ ",MOLECULE,ACTIVE,NULL,2.27e-08); 00145 sp = newspecies("N2 ",MOLECULE,ACTIVE,NULL,2.46e-17); 00146 sp = newspecies("N2+ ",MOLECULE,ACTIVE,NULL,2.22e-17); 00147 sp = newspecies("NO ",MOLECULE,ACTIVE,NULL,5.99e-12); 00148 sp = newspecies("NO+ ",MOLECULE,ACTIVE,NULL,1.31e-11); 00149 sp = newspecies("S2 ",MOLECULE,ACTIVE,NULL,1.21e-11); 00150 sp = newspecies("S2+ ",MOLECULE,ACTIVE,NULL,1.00e-13); 00151 sp = newspecies("OCN ",MOLECULE,ACTIVE,NULL,3.86e-11); 00152 sp = newspecies("OCN+ ",MOLECULE,ACTIVE,NULL,1.65e-22); 00153 sp = newspecies("NH ",MOLECULE,ACTIVE,NULL,1.30e-09); 00154 sp = newspecies("NH+ ",MOLECULE,ACTIVE,NULL,2.20e-10); 00155 sp = newspecies("NH2 ",MOLECULE,ACTIVE,NULL,6.78e-20); 00156 sp = newspecies("NH2+ ",MOLECULE,ACTIVE,NULL,1.71e-16); 00157 sp = newspecies("NH3 ",MOLECULE,ACTIVE,NULL,1.25e-29); 00158 sp = newspecies("NH3+ ",MOLECULE,ACTIVE,NULL,3.03e-23); 00159 sp = newspecies("NH4+ ",MOLECULE,ACTIVE,NULL,1.15e-33); 00160 sp = newspecies("CN ",MOLECULE,ACTIVE,NULL,3.38e-11); 00161 sp = newspecies("CN+ ",MOLECULE,ACTIVE,NULL,5.07e-11); 00162 sp = newspecies("HCN ",MOLECULE,ACTIVE,NULL,1.07e-17); 00163 sp = newspecies("HCN+ ",MOLECULE,ACTIVE,NULL,9.89e-17); 00164 sp = newspecies("HNO ",MOLECULE,ACTIVE,NULL,9.09e-21); 00165 sp = newspecies("HNO+ ",MOLECULE,ACTIVE,NULL,1.91e-18); 00166 sp = newspecies("HS ",MOLECULE,ACTIVE,NULL,1.71e-14); 00167 sp = newspecies("HS+ ",MOLECULE,ACTIVE,NULL,4.83e-13); 00168 sp = newspecies("CS ",MOLECULE,ACTIVE,NULL,6.69e-18); 00169 sp = newspecies("CS+ ",MOLECULE,ACTIVE,NULL,4.12e-11); 00170 sp = newspecies("NO2 ",MOLECULE,ACTIVE,NULL,2.67e-26); 00171 sp = newspecies("NO2+ ",MOLECULE,ACTIVE,NULL,2.41e-23); 00172 sp = newspecies("NS ",MOLECULE,ACTIVE,NULL,3.12e-11); 00173 sp = newspecies("NS+ ",MOLECULE,ACTIVE,NULL,6.81e-13); 00174 sp = newspecies("SO ",MOLECULE,ACTIVE,NULL,4.00e-15); 00175 sp = newspecies("SO+ ",MOLECULE,ACTIVE,NULL,1.20e-07); 00176 sp = newspecies("SiN ",MOLECULE,ACTIVE,NULL,7.03e-12); 00177 sp = newspecies("SiN+ ",MOLECULE,ACTIVE,NULL,7.60e-14); 00178 sp = newspecies("N2O ",MOLECULE,ACTIVE,NULL,1.05e-11); 00179 sp = newspecies("HCS+ ",MOLECULE,ACTIVE,NULL,8.91e-17); 00180 sp = newspecies("OCS ",MOLECULE,ACTIVE,NULL,8.83e-24); 00181 sp = newspecies("OCS+ ",MOLECULE,ACTIVE,NULL,1.72e-21); 00182 sp = newspecies("HNC ",MOLECULE,ACTIVE,NULL,4.00e-15); 00183 sp = newspecies("HCNH+ ",MOLECULE,ACTIVE,NULL,4.00e-15); 00184 sp = newspecies("C2 ",MOLECULE,ACTIVE,NULL,3.71e-09); 00185 sp = newspecies("C2+ ",MOLECULE,ACTIVE,NULL,8.59e-13); 00186 sp = newspecies("CCl ",MOLECULE,ACTIVE,NULL,4.00e-15); 00187 sp = newspecies("ClO ",MOLECULE,ACTIVE,NULL,4.00e-15); 00188 sp = newspecies("HCl+ ",MOLECULE,ACTIVE,NULL,4.00e-15); 00189 sp = newspecies("HCl ",MOLECULE,ACTIVE,NULL,4.00e-15); 00190 sp = newspecies("H2Cl+ ",MOLECULE,ACTIVE,NULL,4.00e-15); 00191 sp = newspecies("CCl+ ",MOLECULE,ACTIVE,NULL,4.00e-15); 00192 sp = newspecies("H2CCl+",MOLECULE,ACTIVE,NULL,0.00e+00); 00193 sp = newspecies("ClO+ ",MOLECULE,ACTIVE,NULL,4.00e-15); 00194 /* fprintf(stderr,"ETC: %d %d\n",gv.lgDustOn, mole.lgGrain_mole_deplete); */ 00195 if(gv.lgDustOn && mole.lgGrain_mole_deplete ) 00196 { 00197 sp = newspecies("COgrn ",MOLECULE,ACTIVE,NULL,1.00e-15); 00198 sp = newspecies("H2Ogrn",MOLECULE,ACTIVE,NULL,0.00e+00); 00199 sp = newspecies("OHgrn ",MOLECULE,ACTIVE,NULL,1.00e-15); 00200 } 00201 sp = newspecies("C2H ",MOLECULE,ACTIVE,NULL,4.00e-15); 00202 sp = newspecies("C2H+ ",MOLECULE,ACTIVE,NULL,4.00e-15); 00203 sp = newspecies("C3 ",MOLECULE,ACTIVE,NULL,1.00e-15); 00204 sp = newspecies("C3+ ",MOLECULE,ACTIVE,NULL,4.00e-15); 00205 sp = newspecies("C3H+ ",MOLECULE,ACTIVE,NULL,4.00e-15); 00206 sp = newspecies("C3H ",MOLECULE,ACTIVE,NULL,4.00e-15); 00207 sp = newspecies("C2H2 ",MOLECULE,ACTIVE,NULL,4.00e-15); 00208 sp = newspecies("C2H2+ ",MOLECULE,ACTIVE,NULL,4.00e-15); 00209 sp = newspecies("C2H3+ ",MOLECULE,ACTIVE,NULL,4.00e-15); 00210 sp = newspecies("N2H+ ",MOLECULE,ACTIVE,NULL,4.00e-15); 00211 /* Add passive species to complete network */ 00212 sp = newspecies("e- ",OTHER,PASSIVE,&(dense.eden_f),0.00e+00); 00213 sp->nElec = -1; sp->mole_mass = (realnum)ELECTRON_MASS; /* Augment properties for this non-molecular species */ 00214 00215 sp = newspecies("Fe ",MOLECULE,PASSIVE,&(dense.xIonDense[ipIRON][0]),0.00e+00); 00216 sp = newspecies("Fe+ ",MOLECULE,PASSIVE,&(dense.xIonDense[ipIRON][1]),0.00e+00); 00217 sp = newspecies("H ",MOLECULE,PASSIVE,&(dense.xIonDense[ipHYDROGEN][0]),0.00e+00); 00218 sp = newspecies("H- ",MOLECULE,PASSIVE,&(hmi.Hmolec[ipMHm]),0.00e+00); 00219 sp = newspecies("H+ ",MOLECULE,PASSIVE,&(dense.xIonDense[ipHYDROGEN][1]),0.00e+00); 00220 sp = newspecies("H2 ",MOLECULE,PASSIVE,&(hmi.Hmolec[ipMH2g]),0.00e+00); 00221 sp = newspecies("H2* ",MOLECULE,PASSIVE,&(hmi.Hmolec[ipMH2s]),0.00e+00); 00222 sp = newspecies("H2+ ",MOLECULE,PASSIVE,&(hmi.Hmolec[ipMH2p]),0.00e+00); 00223 sp = newspecies("H3+ ",MOLECULE,PASSIVE,&(hmi.Hmolec[ipMH3p]),0.00e+00); 00224 sp = newspecies("He ",MOLECULE,PASSIVE,&(dense.xIonDense[ipHELIUM][0]),0.00e+00); 00225 sp = newspecies("He+ ",MOLECULE,PASSIVE,&(dense.xIonDense[ipHELIUM][1]),0.00e+00); 00226 sp = newspecies("Mg ",MOLECULE,PASSIVE,&(dense.xIonDense[ipMAGNESIUM][0]),0.00e+00); 00227 sp = newspecies("Mg+ ",MOLECULE,PASSIVE,&(dense.xIonDense[ipMAGNESIUM][1]),0.00e+00); 00228 00229 /* Create linear list of species and populate it... */ 00230 COmole = (struct molecule **)MALLOC((size_t)mole.num_comole_tot* 00231 sizeof(struct molecule *)); 00232 00233 /* ...first active species */ 00234 i = makeplist(mole_priv.spectab,(void **)COmole, 00235 mole.num_comole_calc,isactive); 00236 ASSERT (i == mole.num_comole_calc); 00237 00238 /* ...then passive species at end of list */ 00239 i = makeplist(mole_priv.spectab,(void **)COmole+mole.num_comole_calc, 00240 mole.num_comole_tot-mole.num_comole_calc,ispassive); 00241 ASSERT (i == mole.num_comole_tot-mole.num_comole_calc); 00242 00243 /* Set molecule indices to order of list just created */ 00244 for(i=0;i<mole.num_comole_tot;i++) 00245 { 00246 COmole[i]->index = i; 00247 } 00248 00249 /* Register the atomic ladders which are active in this network */ 00250 for(i=0;i<mole.num_comole_calc;i++) 00251 { 00252 sp = COmole[i]; 00253 if(sp->active && sp->n_nuclei == 1) 00254 { 00255 if(sp->nElec == 0) 00256 { 00257 element_list[sp->nelem_hevmol]->ipMl = i; 00258 } 00259 else if(sp->nElec == 1) 00260 { 00261 element_list[sp->nelem_hevmol]->ipMlP = i; 00262 } 00263 } 00264 } 00265 00266 /* Create and fill cache for looping on active atomic species */ 00267 chem_element = (struct chem_element_s **) 00268 MALLOC((size_t)(mole.num_elements)*sizeof(struct chem_element_s *)); 00269 i = makeplist(mole_priv.elemtab,(void **)chem_element,mole.num_elements,isCOnet); 00270 ASSERT(i == mole.num_elements); 00271 00272 00273 /* Create other workspace arrays which have sizes given by the species numbers */ 00274 mole.amat = (double **)MALLOC((size_t)mole.num_comole_calc*sizeof(double *)); 00275 mole.amat[0] = (double *)MALLOC((size_t)mole.num_comole_calc* 00276 mole.num_comole_calc*sizeof(double)); 00277 for(i=1;i<mole.num_comole_calc;i++) 00278 { 00279 mole.amat[i] = mole.amat[i-1]+mole.num_comole_calc; 00280 } 00281 mole.b = (double *)MALLOC((size_t)mole.num_comole_calc*sizeof(double)); 00282 mole.c = (double **)MALLOC((size_t)mole.num_comole_tot*sizeof(double *)); 00283 mole.c[0] = (double *)MALLOC((size_t)mole.num_comole_tot* 00284 mole.num_comole_calc*sizeof(double)); 00285 for(i=1;i<mole.num_comole_tot;i++) 00286 { 00287 mole.c[i] = mole.c[i-1]+mole.num_comole_calc; 00288 } 00289 ipiv = (int32 *)MALLOC((size_t)mole.num_comole_calc*sizeof(int32)); 00290 00291 tot_ion = (realnum *)MALLOC((size_t)mole.num_elements*sizeof(realnum)); 00292 00293 return; 00294 00295 } 00296 /*lint +e778 constant expression evaluates to 0 in operation '-' */ 00297 00298 /* Fill element linking structure */ 00299 STATIC void newelement(const char label[], int ipion, int priority) 00300 { 00301 struct chem_element_s *element; 00302 data_u *p; 00303 int exists; 00304 00305 DEBUG_ENTRY("newelement()"); 00306 00307 element = (struct chem_element_s *) MALLOC(sizeof(struct chem_element_s)); 00308 element->ipCl = ipion; 00309 element->ipZ = priority; 00310 ASSERT ((size_t)strlen(label) < sizeof(element->chName)); 00311 strncpy(element->chName,label,sizeof(element->chName)); 00312 element->ipMl = element->ipMlP = -1; /* Chemical network species indices not yet defined */ 00313 p = addentry(element->chName,0,mole_priv.elemtab,&exists); 00314 p->p = (void *) element; 00315 element_list[ipion] = element; 00316 } 00317 /* Parse species string to find constituent atoms, charge etc. */ 00318 #include <string.h> 00319 #include <ctype.h> 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[3], *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(mol->label,' '); 00347 if(s) 00348 *s = '\0'; 00349 00350 /* Insert species into hash table */ 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 /* Copy to private workspace */ 00368 strncpy(mylab,label,7); 00369 00370 /* Trailing modifiers */ 00371 mol->nElec = mol->Excit = 0; 00372 00373 /* Excitation... */ 00374 s = strpbrk(mylab,"*"); 00375 if(s) 00376 { 00377 mol->Excit = 1; 00378 *s = '\0'; 00379 } 00380 00381 /* ...Charge */ 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 /* ...Grain */ 00396 s = strstr(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 /* Now analyse chemical formula */ 00408 nnuc = 0; /* Keep account of number of atoms contained */ 00409 i = 0; 00410 maxel = NULL; 00411 while (mylab[i] != '\0' && mylab[i] != ' ' && mylab[i] != '*') 00412 { 00413 /* Select next element in species, matches regexp [A-Z][a-z]? */ 00414 nel = 0; 00415 thisel[nel++] = mylab[i++]; 00416 if(islower(mylab[i])) 00417 { 00418 thisel[nel++] = mylab[i++]; 00419 } 00420 thisel[nel] = '\0'; 00421 00422 el = findelement(thisel); 00423 if(el == NULL) 00424 { 00425 fprintf(stderr,"Did not recognize element at %s in \"%s \"[%ld]\n",mylab+i,label,i); 00426 cdEXIT( EXIT_FAILURE ); 00427 } 00428 00429 /* Determine 'heaviest' atom in molecular species */ 00430 /* >> chng 06 Sep 25 rjrw: only C, O, Si, N, S, Cl and H count for this */ 00431 if(el->ipZ != 0 && (maxel == NULL || el->ipZ > maxel->ipZ)) 00432 { 00433 maxel = el; 00434 } 00435 00436 if(isdigit(mylab[i])) /* If there is >1 of this atom */ 00437 { 00438 n = 0; 00439 do { 00440 n = 10*n+(long int)(mylab[i]-'0'); 00441 i++; 00442 } while (isdigit(mylab[i])); 00443 } 00444 else 00445 { 00446 n = 1; 00447 } 00448 mol->nElem[el->ipCl] += n; 00449 nnuc += n; 00450 } 00451 00452 mol->n_nuclei = nnuc; 00453 00454 ASSERT(mol->active == 0 || maxel != NULL ); 00455 00456 if(maxel != NULL) 00457 mol->nelem_hevmol = maxel->ipCl; 00458 else 00459 mol->nelem_hevmol = -1; 00460 00461 mol->mole_mass = 0.; 00462 for(i=0;i<LIMELM;++i) 00463 { 00464 mol->mole_mass += 00465 mol->nElem[i]*dense.AtomicWeight[i]*(realnum)ATOMIC_MASS_UNIT; 00466 } 00467 00468 /* If we've found a neutral atomic species active in this network */ 00469 if(mol->n_nuclei == 1 && mol->active && mol->nElec == 0) 00470 { 00471 mole.num_elements++; 00472 mole.lgElem_in_chemistry[mol->nelem_hevmol] = true; 00473 } 00474 } 00475 return mol; 00476 } 00477 STATIC int isactive(data_u *dat) 00478 { 00479 00480 DEBUG_ENTRY("isactive()"); 00481 00482 struct molecule *p = (struct molecule *) (dat->p); 00483 00484 return p->active; 00485 } 00486 STATIC int ispassive(data_u *dat) 00487 { 00488 return !((struct molecule *) dat->p)->active; 00489 } 00490 STATIC int isCOnet(data_u *dat) 00491 { 00492 return (((struct chem_element_s *) dat->p)->ipMl != -1); 00493 } 00494 00495 struct molecule *findspecies(const char buf[]) 00496 { 00497 data_u *p; 00498 00499 DEBUG_ENTRY("findspecies()"); 00500 00501 p = lookup(buf,0,mole_priv.spectab); 00502 00503 if(p) 00504 return (struct molecule *) p->p; 00505 else 00506 return &null_mole; 00507 } 00508 00509 STATIC struct chem_element_s *findelement(const char buf[]) 00510 { 00511 data_u *p; 00512 00513 DEBUG_ENTRY("findelement()"); 00514 00515 p = lookup(buf,0,mole_priv.elemtab); 00516 00517 if(p) 00518 return (struct chem_element_s *) p->p; 00519 else 00520 return NULL; 00521 } 00522 00523 struct COmole_rate_s *CO_findrate_s(const char buf[]) 00524 { 00525 data_u *p; 00526 00527 DEBUG_ENTRY("CO_findrate_s()"); 00528 00529 p = lookup(buf,0,mole_priv.reactab); 00530 00531 if(p) 00532 return (struct COmole_rate_s *) p->p; 00533 else 00534 return NULL; 00535 } 00536 double CO_findrk(const char buf[]) 00537 { 00538 struct COmole_rate_s *rate; 00539 00540 DEBUG_ENTRY("CO_findrk()"); 00541 00542 rate = CO_findrate_s(buf); 00543 00544 if(!rate) 00545 return 0.; 00546 00547 /* check for NaN */ 00548 ASSERT( !isnan( rate->rk ) ); 00549 00550 return rate->rk; 00551 } 00552 double CO_findrate(const char buf[]) 00553 { 00554 struct COmole_rate_s *rate; 00555 double ret; 00556 int i; 00557 00558 DEBUG_ENTRY("CO_findrate()"); 00559 00560 rate = CO_findrate_s(buf); 00561 if(!rate) 00562 { 00563 return 0.; 00564 } 00565 00566 ret = rate->rk; 00567 for(i=0;i<rate->nrates;i++) 00568 ret *= rate->rate_species[i]->hevmol; 00569 return ret; 00570 } 00571 00572 void CO_update_species_cache(void) 00573 { 00574 int i; 00575 00576 DEBUG_ENTRY("CO_update_species_cache()"); 00577 00578 dense.eden_f = (realnum)dense.eden; /* Need floating point version for compatibility with all other values */ 00579 00580 for(i=0;i<mole.num_comole_tot;i++) 00581 { 00582 if(COmole[i]->location) 00583 { 00584 COmole[i]->hevmol = *(COmole[i]->location); 00585 00586 /* check for NaN */ 00587 ASSERT( !isnan( COmole[i]->hevmol ) ); 00588 } 00589 } 00590 } 00591 00592 /* Calculate rate at which molecular network abstracts species */ 00593 00594 /* Need to check rate_species vs reactant behaviour */ 00595 double CO_sink_rate(const char chSpecies[7]) 00596 { 00597 int ipthis, i, n, nt; 00598 double ratev, ratevi; 00599 struct COmole_rate_s *rate; 00600 struct molecule *sp; 00601 00602 DEBUG_ENTRY("CO_sink_rate()"); 00603 00604 sp = findspecies(chSpecies); 00605 ratev = 0; 00606 nt = coreactions.n; 00607 00608 for(n=0;n<nt;n++) 00609 { 00610 rate = coreactions.list[n]; 00611 ipthis = -1; 00612 for(i=0;i<rate->nrates && ipthis == -1;i++) 00613 { 00614 if(rate->rate_species[i] == sp) 00615 { 00616 ipthis = i; 00617 } 00618 } 00619 if(ipthis != -1) { 00620 ratevi = rate->rk; 00621 for(i=0;i<rate->nrates;i++) 00622 { 00623 if(i!=ipthis) 00624 { 00625 ratevi *= rate->rate_species[i]->hevmol; 00626 } 00627 } 00628 ratev += ratevi; 00629 } 00630 } 00631 00632 return ratev; 00633 } 00634 #ifdef PRINTREACT 00635 STATIC void printreact(struct COmole_rate_s *rate) 00636 { 00637 int i; 00638 00639 DEBUG_ENTRY("printreact()"); 00640 00641 for(i=0;i<rate->nreactants;i++) { 00642 fprintf(stderr,"%s,",rate->reactants[i]->label); 00643 } 00644 fprintf(stderr,"=>"); 00645 for(i=0;i<rate->nproducts;i++) { 00646 fprintf(stderr,"%s,",rate->products[i]->label); 00647 } 00648 fprintf(stderr,"\n"); 00649 00650 } 00651 #endif 00652 void CO_update_rks(void) 00653 { 00654 int n, nt; 00655 struct COmole_rate_s *rate; 00656 00657 DEBUG_ENTRY("CO_update_rks()"); 00658 00659 nt = coreactions.n; 00660 00661 for(n=0;n<nt;n++) 00662 { 00663 rate = coreactions.list[n]; 00664 if(rate->fun != NULL) { 00665 rate->rk = rate->a*rate->fun(rate); 00666 } 00667 } 00668 } 00669 00670 double CO_dissoc_rate(const char chSpecies[7]) 00671 { 00672 int ipthis, i, n, nt; 00673 double ratev, ratevi; 00674 struct COmole_rate_s *rate; 00675 struct molecule *sp; 00676 00677 DEBUG_ENTRY("CO_dissoc_rate()"); 00678 00679 sp = findspecies(chSpecies); 00680 ratev = 0; 00681 nt = coreactions.n; 00682 00683 for(n=0;n<nt;n++) 00684 { 00685 rate = coreactions.list[n]; 00686 if(rate->photon == -1) 00687 { 00688 ipthis = 0; 00689 for(i=0;i<rate->nproducts;i++) 00690 { 00691 if(rate->products[i] == sp) 00692 { 00693 ipthis++; 00694 } 00695 } 00696 if(ipthis) 00697 { 00698 ratevi = rate->rk; 00699 for(i=0;i<rate->nrates;i++) 00700 { 00701 ratevi *= rate->rate_species[i]->hevmol; 00702 } 00703 ratev += ipthis*ratevi; 00704 } 00705 } 00706 } 00707 00708 return ratev; 00709 } 00710 double CO_source_rate(const char chSpecies[7]) 00711 { 00712 int ipthis, i, n, nt; 00713 double ratev, ratevi; 00714 struct COmole_rate_s *rate; 00715 struct molecule *sp; 00716 00717 DEBUG_ENTRY("CO_source_rate()"); 00718 00719 sp = findspecies(chSpecies); 00720 ratev = 0; 00721 nt = coreactions.n; 00722 00723 for(n=0;n<nt;n++) 00724 { 00725 rate = coreactions.list[n]; 00726 ipthis = 0; 00727 for(i=0;i<rate->nproducts;i++) 00728 { 00729 if(rate->products[i] == sp) 00730 { 00731 ipthis++; 00732 } 00733 } 00734 if(ipthis) { 00735 ratevi = rate->rk; 00736 for(i=0;i<rate->nrates;i++) 00737 { 00738 ratevi *= rate->rate_species[i]->hevmol; 00739 } 00740 ratev += ipthis*ratevi; 00741 } 00742 } 00743 00744 return ratev; 00745 } 00746 00747 /* Punch all rates involving specified species */ 00748 void CO_punch_mol(FILE *punit, const char chSpecies[], char header[], double depth) 00749 { 00750 int n, nt, i, ipthis; 00751 double ratevi; 00752 struct COmole_rate_s *rate; 00753 struct molecule *sp; 00754 char *s; 00755 00756 DEBUG_ENTRY("CO_punch_mol()"); 00757 00758 s = header; 00759 00760 sp = findspecies(chSpecies); 00761 if(punit == NULL) 00762 { 00763 sprintf (s,"#Depth"); 00764 s += strlen(s); 00765 } 00766 else 00767 { 00768 fprintf (punit,"%.5e",depth); 00769 } 00770 00771 nt = coreactions.n; 00772 for(n=0;n<nt;n++) { 00773 rate = coreactions.list[n]; 00774 ipthis = 0; 00775 for(i=0;i<rate->nrates;i++) 00776 { 00777 if(rate->rate_species[i] == sp) 00778 { 00779 ipthis++; 00780 } 00781 } 00782 for(i=0;i<rate->nproducts;i++) 00783 { 00784 if(rate->products[i] == sp) 00785 { 00786 ipthis++; 00787 } 00788 } 00789 if(ipthis) 00790 { 00791 if(punit == NULL) 00792 { 00793 sprintf(s,"\t%s",rate->label); 00794 s += strlen(s); 00795 } 00796 else 00797 { 00798 ratevi = rate->rk; 00799 for(i=0;i<rate->nrates;i++) 00800 { 00801 ratevi *= rate->rate_species[i]->hevmol; 00802 } 00803 fprintf(punit,"\t%.3e",ratevi); 00804 } 00805 } 00806 } 00807 if(punit == NULL) 00808 { 00809 sprintf(s,"\n"); 00810 } 00811 else 00812 { 00813 fprintf(punit,"\n"); 00814 } 00815 } 00816 00817 void CO_zero(void) 00818 { 00819 long int i; 00820 00821 static bool lgFirstCall = true; 00822 static long int num_comole_calc_MALLOC=-1; 00823 00824 DEBUG_ENTRY("CO_zero()"); 00825 00826 if( lgFirstCall ) 00827 { 00828 /* do one-time malloc of timesc.AgeCOMoleDest */ 00829 timesc.AgeCOMoleDest = (double*)MALLOC( mole.num_comole_calc*sizeof(double) ); 00830 00831 lgFirstCall = false; 00832 num_comole_calc_MALLOC = mole.num_comole_calc; 00833 } 00834 else if( mole.num_comole_calc>num_comole_calc_MALLOC ) 00835 { 00836 /* number of species has increased since last time - this can't happen 00837 * tsuite / programs / comp4 has 95 first time, 98 second time */ 00838 fprintf(ioQQQ,"DISASTER - the number of species in the CO network has increased. This is not allowed.\n"); 00839 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"); 00840 fprintf(ioQQQ,"Sorry.\n"); 00841 cdEXIT(EXIT_FAILURE); 00842 } 00843 00844 for( i=0; i < mole.num_comole_calc; i++ ) 00845 { 00846 timesc.AgeCOMoleDest[i] = 0.; 00847 } 00848 00849 /* largest fraction of atoms in molecules */ 00850 for( i=0; i<mole.num_comole_calc; ++i ) 00851 COmole[i]->xMoleFracMax = 0.; 00852 00853 /* zero out molecular species */ 00854 for( i=0; i < mole.num_comole_tot; i++ ) 00855 { 00856 COmole[i]->hevmol = 0.; 00857 COmole[i]->hevcol = 0.; 00858 } 00859 }