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 /*ParsePowerlawContinuum parse the power law continuum command */ 00004 #include "cddefines.h" 00005 #include "rfield.h" 00006 #include "optimize.h" 00007 #include "input.h" 00008 #include "parse.h" 00009 #include "physconst.h" 00010 00011 void ParsePowerlawContinuum(char *chCard) 00012 { 00013 bool lgEOL; 00014 long int i; 00015 00016 DEBUG_ENTRY( "ParsePowerlawContinuum()" ); 00017 00018 /* power law with cutoff and X-ray continuum */ 00019 strcpy( rfield.chSpType[rfield.nspec], "POWER" ); 00020 00021 /* first parameter is slope of continuum, probably should be negative */ 00022 i = 5; 00023 rfield.slope[rfield.nspec] = FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL); 00024 if( lgEOL ) 00025 { 00026 fprintf( ioQQQ, " There should have been a number on this line. Sorry.\n" ); 00027 cdEXIT(EXIT_FAILURE); 00028 } 00029 00030 if( rfield.slope[rfield.nspec] >= 0. ) 00031 { 00032 fprintf( ioQQQ, " Is the slope of this power law correct?\n" ); 00033 } 00034 00035 /* second optional parameter is high energy cut off */ 00036 rfield.cutoff[rfield.nspec][0] = FFmtRead(chCard,&i,INPUT_LINE_LENGTH, 00037 &lgEOL); 00038 00039 /* no cutoff if eof hit */ 00040 if( lgEOL ) 00041 { 00042 /* no extra parameters at all, so put in extreme cutoffs */ 00043 rfield.cutoff[rfield.nspec][0] = 1e4; 00044 rfield.cutoff[rfield.nspec][1] = 1e-4; 00045 } 00046 else 00047 { 00048 /* first cutoff was present, check for second */ 00049 rfield.cutoff[rfield.nspec][1] = FFmtRead(chCard,&i, 00050 INPUT_LINE_LENGTH,&lgEOL); 00051 if( lgEOL ) 00052 rfield.cutoff[rfield.nspec][1] = 1e-4; 00053 } 00054 00055 /* check that energies were entered in the correct order */ 00056 if( rfield.cutoff[rfield.nspec][1] > rfield.cutoff[rfield.nspec][0] ) 00057 { 00058 fprintf( ioQQQ, " The optional cutoff energies do not appear to be in the correct order. Check Hazy.\n" ); 00059 cdEXIT(EXIT_FAILURE); 00060 } 00061 00062 /* check for keyword KELVIN to interpret cutoff energies as degrees */ 00063 if( nMatch("KELV",chCard) ) 00064 { 00065 /* temps are log if first le 10 */ 00066 if( rfield.cutoff[rfield.nspec][0] <= 10. ) 00067 rfield.cutoff[rfield.nspec][0] = pow(10.,rfield.cutoff[rfield.nspec][0]); 00068 if( rfield.cutoff[rfield.nspec][1] <= 10. ) 00069 rfield.cutoff[rfield.nspec][1] = pow(10.,rfield.cutoff[rfield.nspec][1]); 00070 rfield.cutoff[rfield.nspec][0] /= TE1RYD; 00071 rfield.cutoff[rfield.nspec][1] /= TE1RYD; 00072 } 00073 00074 if( rfield.cutoff[rfield.nspec][0] < 0. || 00075 rfield.cutoff[rfield.nspec][1] < 0. ) 00076 { 00077 fprintf( ioQQQ, " A negative cutoff energy is not physical. Sorry.\n" ); 00078 cdEXIT(EXIT_FAILURE); 00079 } 00080 00081 if( rfield.cutoff[rfield.nspec][1] == 0. && 00082 rfield.slope[rfield.nspec] <= -1. ) 00083 { 00084 fprintf( ioQQQ, " A power-law with this slope, and no low energy cutoff, may have an unphysically large\n brightness temperature in the radio.\n" ); 00085 } 00086 00087 /* vary option */ 00088 if( optimize.lgVarOn ) 00089 { 00090 /* pointer to where to write */ 00091 optimize.nvfpnt[optimize.nparm] = input.nRead; 00092 if( nMatch("ARYB",chCard) ) 00093 { 00094 /* this test is for key "varyb", meaning to vary second parameter 00095 * the cutoff temperature 00096 * this is the number of parameters to feed onto the input line */ 00097 optimize.nvarxt[optimize.nparm] = 2; 00098 /* vary ?? */ 00099 sprintf( optimize.chVarFmt[optimize.nparm], 00100 "POWER LAW %f KELVIN%%f %%f", 00101 rfield.slope[rfield.nspec] ); 00102 /* param is linear scale factor */ 00103 optimize.vparm[0][optimize.nparm] = (realnum)log10(rfield.cutoff[rfield.nspec][0]* 00104 TE1RYD); 00105 optimize.vparm[1][optimize.nparm] = (realnum)log10(rfield.cutoff[rfield.nspec][1]* 00106 TE1RYD); 00107 optimize.vincr[optimize.nparm] = 0.2f; 00108 } 00109 else if( nMatch("ARYC",chCard) ) 00110 { 00111 /* the keyword was "varyc" 00112 * this is the number of parameters to feed onto the input line */ 00113 optimize.nvarxt[optimize.nparm] = 1; 00114 00115 /* vary ?? */ 00116 sprintf( optimize.chVarFmt[optimize.nparm], "POWER LAW%f %f KELVIN %%f )", 00117 rfield.slope[rfield.nspec], 00118 log10(rfield.cutoff[rfield.nspec][0]* TE1RYD) ); 00119 00120 optimize.vparm[0][optimize.nparm] = (realnum)log10(rfield.cutoff[rfield.nspec][1]* 00121 TE1RYD); 00122 optimize.vincr[optimize.nparm] = 0.2f; 00123 } 00124 else 00125 { 00126 /* vary the first parameter only, but still are two more 00127 * this is the number of parameters to feed onto the input line */ 00128 optimize.nvarxt[optimize.nparm] = 3; 00129 strcpy( optimize.chVarFmt[optimize.nparm], 00130 "POWER LAW KELVIN%f %f %f" ); 00131 /* param is slope of the power law */ 00132 optimize.vparm[0][optimize.nparm] = (realnum)rfield.slope[rfield.nspec]; 00133 optimize.vparm[1][optimize.nparm] = (realnum)log10(rfield.cutoff[rfield.nspec][0]* 00134 TE1RYD); 00135 optimize.vparm[2][optimize.nparm] = (realnum)log10(rfield.cutoff[rfield.nspec][1]* 00136 TE1RYD); 00137 optimize.vincr[optimize.nparm] = 0.2f; 00138 } 00139 ++optimize.nparm; 00140 } 00141 00142 /*>>chng 06 nov 10, BUGFIX, nspec was incremented before previous branch 00143 * and so crashed with log10 0 since nspec was beyond set values 00144 * caused fpe domain function error with log 0 00145 * fpe caught by Pavel Abolmasov */ 00146 ++rfield.nspec; 00147 if( rfield.nspec >= LIMSPC ) 00148 { 00149 fprintf( ioQQQ, " Too many continua entered; increase LIMSPC\n" ); 00150 cdEXIT(EXIT_FAILURE); 00151 } 00152 00153 return; 00154 }