Changeset 1740 for branches/newmole/source
- Timestamp:
- 01/12/08 12:50:52 (10 months ago)
- Location:
- branches/newmole/source
- Files:
-
- 23 modified
-
conv_base.cpp (modified) (3 diffs)
-
conv_init_solution.cpp (modified) (2 diffs)
-
conv_ioniz.cpp (modified) (1 diff)
-
conv_itercheck.cpp (modified) (2 diffs)
-
dynamics.cpp (modified) (3 diffs)
-
heat_sum.cpp (modified) (1 diff)
-
iter_end_chk.cpp (modified) (1 diff)
-
iter_startend.cpp (modified) (4 diffs)
-
molcol.cpp (modified) (3 diffs)
-
mole.h (modified) (3 diffs)
-
mole_drive.cpp (modified) (14 diffs)
-
mole_eval_balance.cpp (modified) (2 diffs)
-
mole_newton_step.cpp (modified) (10 diffs)
-
mole_reactions.cpp (modified) (2 diffs)
-
mole_solve.cpp (modified) (2 diffs)
-
mole_species.cpp (modified) (17 diffs)
-
pressure_change.cpp (modified) (1 diff)
-
prt_comment.cpp (modified) (2 diffs)
-
prt_zone.cpp (modified) (1 diff)
-
punch_do.cpp (modified) (2 diffs)
-
radius_increment.cpp (modified) (2 diffs)
-
radius_next.cpp (modified) (5 diffs)
-
rt_tau_reset.cpp (modified) (1 diff)
Legend:
- Unmodified
- Added
- Removed
-
branches/newmole/source/conv_base.cpp
r1739 r1740 160 160 static double SumOTS=0. , OldSumOTS[2]={0.,0.}; 161 161 double save_iso_grnd[NISO][LIMELM]; 162 valarray<realnum> mole_save(mole.num_calc); 162 163 163 164 DEBUG_ENTRY( "ConvBase()" ); … … 229 230 for( i=0; i < mole.num_calc; i++ ) 230 231 { 231 COmole[i]->comole_save = (realnum) COmole[i]->hevmol;232 mole_save[i] = (realnum) mole.list[i]->hevmol; 232 233 } 233 234 … … 926 927 for( i=0; i < mole.num_calc; ++i ) 927 928 { 928 if( fabs( COmole[i]->hevmol-COmole[i]->comole_save)/dense.gas_phase[ipHYDROGEN]-1. >929 if( fabs(mole.list[i]->hevmol-mole_save[i])/dense.gas_phase[ipHYDROGEN]-1. > 929 930 conv.EdenErrorAllowed/2.) 930 931 { 931 932 conv.lgConvIoniz = false; 932 sprintf( conv.chConvIoniz, "ch %-4.4s", COmole[i]->label );933 conv.BadConvIoniz[0] = COmole[i]->comole_save/dense.gas_phase[ipHYDROGEN];934 conv.BadConvIoniz[1] = COmole[i]->hevmol/dense.gas_phase[ipHYDROGEN];933 sprintf( conv.chConvIoniz, "ch %-4.4s",mole.list[i]->label ); 934 conv.BadConvIoniz[0] = mole_save[i]/dense.gas_phase[ipHYDROGEN]; 935 conv.BadConvIoniz[1] = mole.list[i]->hevmol/dense.gas_phase[ipHYDROGEN]; 935 936 } 936 937 } -
branches/newmole/source/conv_init_solution.cpp
r1739 r1740 132 132 for (nelem=ipHYDROGEN;nelem<LIMELM;nelem++) 133 133 { 134 if( dense.lgElmtOn[nelem] && COmole[i]->location == NULL )134 if( dense.lgElmtOn[nelem] && mole.list[i]->location == NULL ) 135 135 { 136 realnum dens_elemsp = (realnum) COmole[i]->hevmol*COmole[i]->nElem[nelem];136 realnum dens_elemsp = (realnum) mole.list[i]->hevmol*mole.list[i]->nElem[nelem]; 137 137 if ( FracMoleMax*dense.gas_phase[nelem] < dens_elemsp ) 138 138 { … … 148 148 for(i=0;i<mole.num_calc;++i) 149 149 { 150 OxyInGrains += (1 - COmole[i]->lgGas_Phase)*COmole[i]->hevmol*COmole[i]->nElem[ipOXYGEN];150 OxyInGrains += (1 - mole.list[i]->lgGas_Phase)*mole.list[i]->hevmol*mole.list[i]->nElem[ipOXYGEN]; 151 151 } 152 152 /* this is now fraction of O in ices */ -
branches/newmole/source/conv_ioniz.cpp
r1739 r1740 108 108 dense.eden, 109 109 gv.TotalEden, 110 co.comole_eden,110 mole.elec, 111 111 dense.EdenTrue); 112 112 } -
branches/newmole/source/conv_itercheck.cpp
r1739 r1740 152 152 { 153 153 double differ; 154 if( COmole[i]->n_nuclei == 1)154 if(mole.list[i]->n_nuclei == 1) 155 155 continue; 156 156 157 157 /* was the species abundance and changing? */ 158 differ = (double)fabs( COmole[i]->hevcol_old-COmole[i]->hevcol) /159 (double)SDIV( COmole[i]->hevcol);160 if( ( COmole[i]->hevcol/colden.colden[ipCOL_HTOT] > 1e-5) &&158 differ = (double)fabs(mole.list[i]->hevcol_old-mole.list[i]->hevcol) / 159 (double)SDIV(mole.list[i]->hevcol); 160 if( (mole.list[i]->hevcol/colden.colden[ipCOL_HTOT] > 1e-5) && 161 161 (differ > conv.autocv) ) 162 162 { … … 168 168 strcpy( conv.chNotConverged, "CO mol col" ); 169 169 /*fprintf(ioQQQ,"debugggreset\t CO mole %li %li %.2e %.2e\n", 170 i,iteration, COmole[i]->hevcol_old,COmole[i]->hevcol);*/170 i,iteration,mole.list[i]->hevcol_old,mole.list[i]->hevcol);*/ 171 171 172 172 if( punch.lgPunConv ) 173 173 { 174 174 fprintf( punch.ipPunConv, "CO mol col, old:%.3e new:%.3e\n" , 175 COmole[i]->hevcol_old ,176 COmole[i]->hevcol );175 mole.list[i]->hevcol_old , 176 mole.list[i]->hevcol ); 177 177 } 178 178 } -
branches/newmole/source/dynamics.cpp
r1739 r1740 740 740 for(i=0;i<mole.num_calc;i++) 741 741 { 742 sumh += (realnum) COmole[i]->hevmol*COmole[i]->nElem[nelem];742 sumh += (realnum) mole.list[i]->hevmol*mole.list[i]->nElem[nelem]; 743 743 } 744 744 fprintf(ioQQQ, … … 951 951 frac_next; 952 952 /* Externally represented species are already counted in UpstreamElem */ 953 if( COmole[mol]->location == NULL)953 if(mole.list[mol]->location == NULL) 954 954 { 955 955 for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem ) 956 956 { 957 UpstreamElem[nelem] += Upstream_molecules[mol]* COmole[mol]->nElem[nelem];957 UpstreamElem[nelem] += Upstream_molecules[mol]*mole.list[mol]->nElem[nelem]; 958 958 } 959 959 } … … 983 983 { 984 984 Upstream_molecules[mol] = Old_molecules[ipUpstream][mol]/Old_density[ipUpstream]; 985 if( COmole[mol]->location == NULL)985 if(mole.list[mol]->location == NULL) 986 986 { 987 987 for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem ) 988 988 { 989 UpstreamElem[nelem] += Upstream_molecules[mol]* COmole[mol]->nElem[nelem];989 UpstreamElem[nelem] += Upstream_molecules[mol]*mole.list[mol]->nElem[nelem]; 990 990 } 991 991 } -
branches/newmole/source/heat_sum.cpp
r1739 r1740 116 116 for( i=0; i < mole.num_calc; i++ ) 117 117 { 118 if( COmole[i]->location == NULL)119 AtomicCollidDensity += COmole[i]->hevmol;118 if(mole.list[i]->location == NULL) 119 AtomicCollidDensity += mole.list[i]->hevmol; 120 120 } 121 121 -
branches/newmole/source/iter_end_chk.cpp
r1739 r1740 51 51 { 52 52 /* define the abundance of oxygen frozen out on grain surfaces */ 53 oxy_in_grains += (1 - COmole[i]->lgGas_Phase)*COmole[i]->hevmol*COmole[i]->nElem[ipOXYGEN];53 oxy_in_grains += (1 - mole.list[i]->lgGas_Phase)*mole.list[i]->hevmol*mole.list[i]->nElem[ipOXYGEN]; 54 54 } 55 55 /*fprintf(ioQQQ, "DEBUG oxy in grains %.2e %e %e\n", -
branches/newmole/source/iter_startend.cpp
r1739 r1740 298 298 for( i=0; i<mole.num_calc; i++) 299 299 { 300 COmole[i]->hevmol_save = (realnum) COmole[i]->hevmol;300 mole.list[i]->hevmol_save = (realnum) mole.list[i]->hevmol; 301 301 } 302 302 ortho_save = h2.ortho_density; … … 433 433 for( i=0; i < mole.num_calc; i++ ) 434 434 { 435 COmole[i]->hevcol = 0.;435 mole.list[i]->hevcol = 0.; 436 436 } 437 437 for( i=0; i < 5; i++ ) … … 890 890 for( i=0; i<mole.num_calc; i++) 891 891 { 892 COmole[i]->hevmol = COmole[i]->hevmol_save;892 mole.list[i]->hevmol = mole.list[i]->hevmol_save; 893 893 } 894 894 hmi.H2_total = (realnum) (findspecies("H2")->hevmol + findspecies("H2*")->hevmol); … … 908 908 { 909 909 /* column densities */ 910 COmole[i]->hevcol = 0.;910 mole.list[i]->hevcol = 0.; 911 911 /* largest fraction of atoms in molecules */ 912 COmole[i]->xMoleFracMax = 0.;913 COmole[i]->xMoleFracMaxElem = -1;912 mole.list[i]->xMoleFracMax = 0.; 913 mole.list[i]->xMoleFracMaxElem = -1; 914 914 } 915 915 -
branches/newmole/source/molcol.cpp
r1739 r1740 57 57 for( i=0; i < mole.num_calc; i++ ) 58 58 { 59 COmole[i]->hevcol = 0.;59 mole.list[i]->hevcol = 0.; 60 60 } 61 61 } … … 66 66 for( i=0; i < mole.num_calc; i++ ) 67 67 { 68 COmole[i]->hevcol += (realnum) COmole[i]->hevmol*radius.drad_x_fillfac;68 mole.list[i]->hevcol += (realnum) mole.list[i]->hevmol*radius.drad_x_fillfac; 69 69 } 70 70 } … … 79 79 for( i=0; i < mole.num_calc; i++ ) 80 80 { 81 if( COmole[i]->location == NULL)81 if(mole.list[i]->location == NULL) 82 82 { 83 fprintf( ioMEAN, " %-6.6s:", COmole[i]->label );84 fprintf( ioMEAN, "%7.3f",log10(MAX2(SMALLFLOAT, COmole[i]->hevcol )));83 fprintf( ioMEAN, " %-6.6s:", mole.list[i]->label ); 84 fprintf( ioMEAN, "%7.3f",log10(MAX2(SMALLFLOAT,mole.list[i]->hevcol ))); 85 85 j++; 86 86 if( j == 8 ) -
branches/newmole/source/mole.h
r1739 r1740 61 61 can be important to the ionization balance */ 62 62 63 struct molecule; 64 63 65 class Mole { 64 66 … … 121 123 *snk; 122 124 realnum ***xMoleChTrRate;/***[LIMELM][LIMELM+1][LIMELM+1];*/ 125 126 127 struct molecule **list; 123 128 }; 124 129 … … 149 154 realnum hev_reinit; /* this saves solution in first zone of previous iteration */ 150 155 realnum hevcol_old; /** column density in previous iteration */ 151 realnum comole_save;152 156 realnum hevmol_save; 153 157 enum mole_state state; 154 158 int index, groupnum; 155 159 }; 156 157 extern struct molecule **COmole;158 160 159 161 extern struct molecule *findspecies(const char buf[]); -
branches/newmole/source/mole_drive.cpp
r1687 r1740 42 42 STATIC void mole_h_rate_diagnostics(void); 43 43 44 struct molecule **COmole;45 46 47 44 /*mole_drive main driver for heavy molecular equilibrium routines */ 48 45 void mole_drive(void) … … 132 129 for( i=0; i < mole.num_calc; i++ ) 133 130 { 134 COmole[i]->hevmol = COmole[i]->hev_reinit*ratio;131 mole.list[i]->hevmol = mole.list[i]->hev_reinit*ratio; 135 132 } 136 133 … … 150 147 for( i=0; i < mole.num_calc; i++ ) 151 148 { 152 COmole[i]->hev_reinit = (realnum) COmole[i]->hevmol;149 mole.list[i]->hev_reinit = (realnum) mole.list[i]->hevmol; 153 150 } 154 151 … … 177 174 for( i=0; i<mole.num_calc; ++i ) 178 175 { 179 COmole[i]->hevmol = 0.;176 mole.list[i]->hevmol = 0.; 180 177 } 181 178 /* heating due to CO photodissociation */ … … 508 505 hmi.HeatH2Dexc_ELWERT = ((mole_findrate("H2*,H2=>H2,H2") + mole_findrate("H2*,H=>H2,H")) 509 506 - (mole_findrate("H2,H2=>H2*,H2") + mole_findrate("H2,H=>H2*,H"))) * 4.17e-12 * f3 + 510 f4 * ( COmole[element_list[ipHYDROGEN]->ipMl[0]]->hevmol/dense.gas_phase[ipHYDROGEN]) +507 f4 * (mole.list[element_list[ipHYDROGEN]->ipMl[0]]->hevmol/dense.gas_phase[ipHYDROGEN]) + 511 508 f5 * secondaries.x12tot * EN1EV * hmi.H2_total; 512 509 … … 788 785 { 789 786 if (nelem != ipHYDROGEN && nelem != ipHELIUM) 790 sum_first_ions += COmole[element_list[nelem]->ipMl[1]]->hevmol;787 sum_first_ions += mole.list[element_list[nelem]->ipMl[1]]->hevmol; 791 788 } 792 789 … … 798 795 destroy_coll_heavies = mole_findrk("H-,Li+=>H,Li")*sum_first_ions; 799 796 destroy_coll_electrons = mole_findrk("H-,e-=>H-,e-,e-")*findspecies("e-")->hevmol; 800 destroy_Hattach = COmole[element_list[ipHYDROGEN]->ipMl[0]]->hevmol*(mole_findrk("H,H-=>H2,e-")+mole_findrk("H,H-=>H2*,e-"));801 destroy_fhneut = mole_findrk("H-,H+=>H,H")* COmole[element_list[ipHYDROGEN]->ipMl[1]]->hevmol;797 destroy_Hattach = mole.list[element_list[ipHYDROGEN]->ipMl[0]]->hevmol*(mole_findrk("H,H-=>H2,e-")+mole_findrk("H,H-=>H2*,e-")); 798 destroy_fhneut = mole_findrk("H-,H+=>H,H")*mole.list[element_list[ipHYDROGEN]->ipMl[1]]->hevmol; 802 799 803 800 destsum = destroy_photo + destroy_coll_heavies + destroy_coll_electrons + … … 835 832 for( i=0; i < mole.num_calc; i++ ) 836 833 { 837 fprintf( ioQQQ, "\t%-4.4s\t%.2e", COmole[i]->label, COmole[i]->hevmol );834 fprintf( ioQQQ, "\t%-4.4s\t%.2e", mole.list[i]->label, mole.list[i]->hevmol ); 838 835 } 839 836 fprintf( ioQQQ, " \n" ); … … 847 844 for( i=0; i < mole.num_calc; i++ ) 848 845 { 849 fprintf( ioQQQ, " %-4.4s:%.2e", COmole[i]->label, COmole[i]->hevmol );846 fprintf( ioQQQ, " %-4.4s:%.2e", mole.list[i]->label, mole.list[i]->hevmol ); 850 847 } 851 848 fprintf( ioQQQ, " \n" ); … … 904 901 { 905 902 if (nelem != ipHYDROGEN && nelem != ipHELIUM) 906 sum_first_ions += COmole[element_list[nelem]->ipMl[1]]->hevmol;903 sum_first_ions += mole.list[element_list[nelem]->ipMl[1]]->hevmol; 907 904 } 908 905 … … 917 914 hmi.HMinus_photo_rate, mole_findrk("H-,Li+=>H,Li")*sum_first_ions, 918 915 mole_findrk("H-,e-=>H-,e-,e-")*findspecies("e-")->hevmol, 919 COmole[element_list[ipHYDROGEN]->ipMl[0]]->hevmol*(mole_findrk("H,H-=>H2,e-")+mole_findrk("H,H-=>H2*,e-")), iso.gamnc[ipH_LIKE][ipHYDROGEN][ipH1s],920 COmole[element_list[ipHYDROGEN]->ipMl[1]]->hevmol*mole_findrk("H-,H+=>H,H") );916 mole.list[element_list[ipHYDROGEN]->ipMl[0]]->hevmol*(mole_findrk("H,H-=>H2,e-")+mole_findrk("H,H-=>H2*,e-")), iso.gamnc[ipH_LIKE][ipHYDROGEN][ipH1s], 917 mole.list[element_list[ipHYDROGEN]->ipMl[1]]->hevmol*mole_findrk("H-,H+=>H,H") ); 921 918 fprintf( ioQQQ, " H- heating:%10.3e Ind cooling%10.2e Spon cooling%10.2e\n", 922 919 hmi.hmihet, hmi.HMinus_induc_rec_cooling, hmi.hmicol ); … … 993 990 /*@-redef@*/ 994 991 /* often the H- route is the most efficient formation mechanism for H2, 995 * will be through rate called COmole[element_list[ipHYDROGEN]->ipMl]->hevmol*hmi.assoc_detach992 * will be through rate called mole.list[element_list[ipHYDROGEN]->ipMl]->hevmol*hmi.assoc_detach 996 993 * this debug print statement is to trace h2 oscillations */ 997 994 enum {DEBUG_LOC=false}; … … 1011 1008 mole_findrk("H,e-=>H-,PHOTON")*findspecies("e-")->hevmol, 1012 1009 hmi.HMinus_induc_rec_rate, 1013 COmole[element_list[ipHYDROGEN]->ipMl[0]]->hevmol,1014 COmole[element_list[ipHYDROGEN]->ipMl[1]]->hevmol,1010 mole.list[element_list[ipHYDROGEN]->ipMl[0]]->hevmol, 1011 mole.list[element_list[ipHYDROGEN]->ipMl[1]]->hevmol, 1015 1012 findspecies("H-")->hevmol, 1016 1013 hmi.H2_total, … … 1056 1053 fprintf(ioQQQ,"Temperature %g\n",phycon.te); 1057 1054 fprintf(ioQQQ," Net mol ion rate [%g %g] %g\n",mole.source[ipHYDROGEN][1],mole.sink[ipHYDROGEN][1], 1058 mole.source[ipHYDROGEN][1]-mole.sink[ipHYDROGEN][1]* COmole[element_list[ipHYDROGEN]->ipMl[1]]->hevmol);1055 mole.source[ipHYDROGEN][1]-mole.sink[ipHYDROGEN][1]*mole.list[element_list[ipHYDROGEN]->ipMl[1]]->hevmol); 1059 1056 } 1060 1057 } -
branches/newmole/source/mole_eval_balance.cpp
r1712 r1740 170 170 for (j=0;j<mole.num_calc;j++) 171 171 { 172 test += c[i][j]* COmole[j]->nElem[nelem];173 tot += fabs(c[i][j])* COmole[j]->nElem[nelem];172 test += c[i][j]*mole.list[j]->nElem[nelem]; 173 tot += fabs(c[i][j])*mole.list[j]->nElem[nelem]; 174 174 } 175 175 if (fabs(test) > 1e-10*tot) 176 176 { 177 177 fprintf(stdout,"Network conservation error %s %s %g %g %g %g\n", 178 COmole[element_list[nelem]->ipMl[0]]->label,COmole[i]->label,179 test,test/tot, COmole[element_list[nelem]->ipMl[0]]->hevmol,COmole[element_list[nelem]->ipMl[1]]->hevmol);178 mole.list[element_list[nelem]->ipMl[0]]->label,mole.list[i]->label, 179 test,test/tot,mole.list[element_list[nelem]->ipMl[0]]->hevmol,mole.list[element_list[nelem]->ipMl[1]]->hevmol); 180 180 //fprintf(stdout,"Problem at %s\n",rate->label); 181 181 exit(-1); … … 188 188 for (i=0;i<num_total;i++) 189 189 { 190 b[i] = mole.src[i]-mole.snk[i]* COmole[i]->hevmol;190 b[i] = mole.src[i]-mole.snk[i]*mole.list[i]->hevmol; 191 191 } 192 192 -
branches/newmole/source/mole_newton_step.cpp
r1714 r1740 118 118 for(i=0;i<mole.num_total;i++) 119 119 { 120 calcv[i] = COmole[i]->hevmol;120 calcv[i] = mole.list[i]->hevmol; 121 121 } 122 122 … … 210 210 for(i=0;i<mole.num_calc;i++) 211 211 { 212 ASSERT( COmole[i]->index == i);212 ASSERT(mole.list[i]->index == i); 213 213 sum = 0.; 214 214 for (ion=0;ion<N_MOLE_ION;ion++) … … 276 276 for(i=0;i<mole.num_calc;i++) 277 277 { 278 c[i][element_list[nelem]->ipMl[0]] = COmole[i]->nElem[nelem]*scale;278 c[i][element_list[nelem]->ipMl[0]] = mole.list[i]->nElem[nelem]*scale; 279 279 } 280 280 } … … 316 316 for( i=0; i < mole.num_calc; i++ ) 317 317 { 318 fprintf( ioQQQ, " %-4.4s", COmole[i]->label );318 fprintf( ioQQQ, " %-4.4s", mole.list[i]->label ); 319 319 } 320 320 fprintf( ioQQQ, " \n" ); … … 322 322 fprintf(ioQQQ," MOLE old abundances\t%.2f",fnzone); 323 323 for( i=0; i<mole.num_calc; i++ ) 324 fprintf(ioQQQ,"\t%.2e", COmole[i]->hevmol );324 fprintf(ioQQQ,"\t%.2e", mole.list[i]->hevmol ); 325 325 fprintf(ioQQQ,"\n" ); 326 326 327 327 for( i=0; i < mole.num_calc; i++ ) 328 328 {