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 /*ParseAgn parse parameters for the AGN continuum shape command */ 00004 #include "cddefines.h" 00005 #include "rfield.h" 00006 #include "parse.h" 00007 #include "physconst.h" 00008 00009 void ParseAgn(char *chCard) 00010 { 00011 bool lgEOL; 00012 long int i; 00013 double BigBump, 00014 Ratio, 00015 XRays, 00016 xnu; 00017 00018 DEBUG_ENTRY( "ParseAgn()" ); 00019 00020 /* this radiation field will be something like an AGN */ 00021 strcpy( rfield.chSpType[rfield.nspec], "AGN " ); 00022 00023 /* there were no numbers on the line - this could be the Kirk option, 00024 * to use the numbers for the continuum in the atlas paper */ 00025 if( nMatch("KIRK",chCard) ) 00026 { 00027 /* million degree cutoff, but in rydbergs */ 00028 rfield.slope[rfield.nspec] = 1e6 / TE1RYD; 00029 00030 /* cutoff is second parameter is really alpha ox */ 00031 rfield.cutoff[rfield.nspec][0] = -1.40; 00032 00033 /* bb slope is third parameter */ 00034 rfield.cutoff[rfield.nspec][1] = -0.50; 00035 00036 /* slope of X-Ray component is last parameter */ 00037 rfield.cutoff[rfield.nspec][2] = -1.0; 00038 } 00039 else 00040 { 00041 /* first parameter is temperature of big bump 00042 * second parameter is alpha ox */ 00043 i = 5; 00044 /* slope is first parameter is really temperature of bump */ 00045 rfield.slope[rfield.nspec] = FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL); 00046 if( lgEOL ) 00047 { 00048 00049 fprintf( ioQQQ, " The big bump temperature should have been on this line. Sorry.\n" ); 00050 cdEXIT(EXIT_FAILURE); 00051 } 00052 00053 if( rfield.slope[rfield.nspec] <= 0. ) 00054 { 00055 fprintf( ioQQQ, " Non positive temperature not allowed. Sorry.\n" ); 00056 cdEXIT(EXIT_FAILURE); 00057 } 00058 00059 /* temps are log if first le 10 */ 00060 if( rfield.slope[rfield.nspec] <= 10. ) 00061 rfield.slope[rfield.nspec] = 00062 pow(10.,rfield.slope[rfield.nspec]); 00063 00064 /* want cutoff in ryd not kelvin */ 00065 rfield.slope[rfield.nspec] /= TE1RYD; 00066 00067 /* cutoff is second parameter is really alpha ox */ 00068 rfield.cutoff[rfield.nspec][0] = FFmtRead(chCard,&i,INPUT_LINE_LENGTH, &lgEOL); 00069 if( lgEOL ) 00070 { 00071 fprintf( ioQQQ, " alpha ox should have been on this line. Sorry.\n" ); 00072 cdEXIT(EXIT_FAILURE); 00073 } 00074 00075 if( rfield.cutoff[rfield.nspec][0] > 3. || 00076 rfield.cutoff[rfield.nspec][0] < -3. ) 00077 { 00078 fprintf( ioQQQ, " An alpha ox of%10.2e looks funny to me. Check Hazy to make sure its ok.\n", 00079 rfield.cutoff[rfield.nspec][0] ); 00080 } 00081 00082 if( rfield.cutoff[rfield.nspec][0] >= 0. ) 00083 { 00084 fprintf( ioQQQ, " The sign of alpha ox is almost certainly incorrect. Check Hazy.\n" ); 00085 } 00086 00087 /* bb slope is third parameter */ 00088 rfield.cutoff[rfield.nspec][1] = FFmtRead(chCard,&i,INPUT_LINE_LENGTH, &lgEOL); 00089 if( lgEOL ) 00090 rfield.cutoff[rfield.nspec][1] = -0.5f; 00091 00092 /* slope of X-Ray component is last parameter */ 00093 rfield.cutoff[rfield.nspec][2] = FFmtRead(chCard,&i,INPUT_LINE_LENGTH, &lgEOL); 00094 if( lgEOL ) 00095 rfield.cutoff[rfield.nspec][2] = -1.0f; 00096 } 00097 00098 /* 403.3 is ratio of energies where alpha ox defined, 00099 * assumed to be 2500A and 2keV */ 00100 Ratio = pow(403.3,rfield.cutoff[rfield.nspec][0] - 1.); 00101 00102 /* following code must be kept parallel with ffun1 */ 00103 xnu = 0.3645; 00104 BigBump = pow(xnu,-1. + rfield.cutoff[rfield.nspec][1])* 00105 sexp(xnu/rfield.slope[rfield.nspec]); 00106 xnu = 147.; 00107 00108 /* XRays = xnu**(-2.) */ 00109 XRays = pow(xnu,rfield.cutoff[rfield.nspec][2] - 1.); 00110 if( BigBump <= 0. ) 00111 { 00112 fprintf( ioQQQ, " Big Bump had zero flux at .3645 Ryd.\n" ); 00113 cdEXIT(EXIT_FAILURE); 00114 } 00115 rfield.cutoff[rfield.nspec][0] = (Ratio/(XRays/BigBump)); 00116 00117 /* lastly increment number of spectra and check that we still have room in the vector */ 00118 ++rfield.nspec; 00119 if( rfield.nspec >= LIMSPC ) 00120 { 00121 fprintf( ioQQQ, " Too many continua entered; increase LIMSPC\n" ); 00122 cdEXIT(EXIT_FAILURE); 00123 } 00124 return; 00125 }