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 /*ParseFluc parse the fluctuations command, which affects either density or abundances */ 00004 #include "cddefines.h" 00005 #include "dense.h" 00006 #include "parse.h" 00007 00008 void ParseFluc(char *chCard ) 00009 { 00010 bool lgEOL; 00011 long int i; 00012 double flmax, 00013 flmin, 00014 period, 00015 temp; 00016 00017 DEBUG_ENTRY( "ParseFluc()" ); 00018 00019 /* rapid density fluctuations 00020 * first parameter is log of period, 2 is log den max, 3 log Nmin */ 00021 if( nMatch("ABUN",chCard) ) 00022 { 00023 /* abundances varied, not density */ 00024 dense.lgDenFlucOn = false; 00025 } 00026 else 00027 { 00028 /* density is varied */ 00029 dense.lgDenFlucOn = true; 00030 } 00031 00032 /* optional keyword COLUMN makes sin over column density rather than radius */ 00033 if( nMatch("COLU",chCard) ) 00034 { 00035 /* found key, not fluc over radius, over col den instead */ 00036 dense.lgDenFlucRadius = false; 00037 } 00038 else 00039 { 00040 /* no key, use default of radius */ 00041 dense.lgDenFlucRadius = true; 00042 } 00043 00044 i = 5; 00045 /* 1st number log of period in centimeters */ 00046 period = pow(10.,FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL)); 00047 dense.flong = (realnum)(6.2831853/period); 00048 temp = FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL); 00049 00050 /* check size of density - will we crash? */ 00051 if( temp > log10(FLT_MAX) || temp < log10(FLT_MIN) ) 00052 { 00053 fprintf(ioQQQ, 00054 " DISASTER - the log of the entered max hydrogen density is %.3f - too extreme for this processor.\n", 00055 temp); 00056 if( temp > 0. ) 00057 fprintf(ioQQQ, 00058 " DISASTER - the log of the largest hydrogen density this processor can do is %.3f.\n", 00059 log10(FLT_MAX) ); 00060 else 00061 fprintf(ioQQQ, 00062 " DISASTER - the log of the smallest hydrogen density this processor can do is %.3f.\n", 00063 log10(FLT_MIN) ); 00064 fprintf(ioQQQ," Sorry.\n" ); 00065 cdEXIT(EXIT_FAILURE); 00066 } 00067 00068 /* 2nd number log of max hydrogen density */ 00069 flmax = pow(10.,temp); 00070 00071 temp = FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL); 00072 00073 /* check size of density - will we crash? */ 00074 if( temp > log10(FLT_MAX) || temp < log10(FLT_MIN) ) 00075 { 00076 fprintf(ioQQQ, 00077 " DISASTER - the log of the entered min hydrogen density is %.3f - too extreme for this processor.\n", 00078 temp); 00079 if( temp > 0. ) 00080 fprintf(ioQQQ, 00081 " DISASTER - the log of the largest hydrogen density this processor can do is %.3f.\n", 00082 log10(FLT_MAX) ); 00083 else 00084 fprintf(ioQQQ, 00085 " DISASTER - the log of the smallest hydrogen density this processor can do is %.3f.\n", 00086 log10(FLT_MIN) ); 00087 fprintf(ioQQQ," Sorry.\n" ); 00088 cdEXIT(EXIT_FAILURE); 00089 } 00090 00091 /* 3rd number log of min hydrogen density */ 00092 flmin = pow(10.,temp); 00093 00094 if( flmax/flmin > 100. ) 00095 { 00096 fprintf( ioQQQ, "This range of density probably will not work.\n" ); 00097 } 00098 if( flmax > 1e15 ) 00099 { 00100 fprintf( ioQQQ, "These parameters look funny to me. Please check Hazy.\n" ); 00101 } 00102 if( lgEOL || (flmin > flmax) ) 00103 { 00104 fprintf( ioQQQ, "There MUST be three numbers on this line.\n" ); 00105 fprintf( ioQQQ, "These must be the period(cm), max, min densities, all logs, in that order.\n" ); 00106 if( flmin > flmax ) 00107 fprintf( ioQQQ, "The max density must be greater or equal than the min density.\n" ); 00108 cdEXIT(EXIT_FAILURE); 00109 } 00110 00111 /* this is optional phase shift for the command */ 00112 dense.flcPhase = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL); 00113 00114 /* FacAbunSav = (cfirst * COS( depth*flong+flcPhase ) + csecnd) */ 00115 dense.cfirst = (realnum)((flmax - flmin)/2.); 00116 dense.csecnd = (realnum)((flmax + flmin)/2.); 00117 /* these will be added together with the first mult by sin - which goes to 00118 * -1 - must not have a negative density */ 00119 ASSERT( dense.cfirst < dense.csecnd ); 00120 /* >>chng 96 jul 13 moved depset to SetAbundances fac 00121 * if( lgDenFlucOn ) then 00122 * this is a pressure law 00123 * chCPres = 'SINE' 00124 * else 00125 * this is the metallicity of the gas 00126 * do i=3,limelm 00127 * depset(i) = flmax 00128 * end do 00129 * endif 00130 * 00131 * now get density if this is density option (not abundances) */ 00132 if( dense.lgDenFlucOn ) 00133 { 00134 strcpy( dense.chDenseLaw, "SINE" ); 00135 00136 if( dense.gas_phase[ipHYDROGEN] > 0. ) 00137 { 00138 fprintf( ioQQQ, " PROBLEM DISASTER More than one density command was entered.\n" ); 00139 cdEXIT(EXIT_FAILURE); 00140 } 00141 00142 /* depth is zero for first zone */ 00143 dense.gas_phase[ipHYDROGEN] = dense.cfirst*(realnum)cos(dense.flcPhase) + dense.csecnd; 00144 00145 if( dense.gas_phase[ipHYDROGEN] <= 0. ) 00146 { 00147 fprintf( ioQQQ, " PROBLEM DISASTER Hydrogen density must be > 0.\n" ); 00148 cdEXIT(EXIT_FAILURE); 00149 } 00150 } 00151 return; 00152 }