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 /*HeLikeError fills uncertainty arrays */ 00004 #include "cddefines.h" 00005 #include "iso.h" 00006 00007 /* This routine handles errors when that option is turned on (via the command 00008 * "atom he-like error generation" */ 00009 void iso_put_error(long int ipISO, 00010 long int nelem, 00011 long int ipHi, 00012 long int ipLo, 00013 long int whichData, 00014 realnum errorToPut) 00015 { 00016 00017 DEBUG_ENTRY( "iso_put_error()" ); 00018 00019 if( iso.lgRandErrGen[ipISO] ) 00020 { 00021 /* whichData is either IPRAD, IPCOLLIS, or IPENERGY */ 00022 ASSERT( whichData <= 2 ); 00023 ASSERT( ipISO <= NISO ); 00024 ASSERT( nelem < LIMELM ); 00025 ASSERT( ipHi <= iso.numLevels_max[ipISO][nelem] ); 00026 ASSERT( ipLo <= iso.numLevels_max[ipISO][nelem] ); 00027 ASSERT( errorToPut >= 0. ); 00028 00029 iso.Error[ipISO][nelem][ipHi][ipLo][whichData] = errorToPut; 00030 } 00031 return; 00032 } 00033 00034 void iso_error_generation( long ipISO, long nelem ) 00035 { 00036 00037 long ipHi, ipLo, typeOfRate; 00038 realnum ErrorToPut; 00039 00040 DEBUG_ENTRY( "iso_error_generation()" ); 00041 00042 /* put recombination errors into iso.Error[ipISO] array */ 00043 for( ipHi=0; ipHi<iso.numLevels_max[ipISO][nelem]; ipHi++ ) 00044 { 00045 if( N_(ipHi)<=5 ) 00046 { 00047 ErrorToPut = 0.02f; 00048 } 00049 else if( N_(ipHi)>iso.n_HighestResolved_max[ipISO][nelem] ) 00050 { 00051 ErrorToPut = 0.02f; 00052 } 00053 else if( L_(ipHi)>=3 ) 00054 { 00055 ErrorToPut = 0.005f; 00056 } 00057 else if( L_(ipHi)==2 ) 00058 { 00059 ErrorToPut = 0.01f; 00060 } 00061 else if( L_(ipHi)==1 ) 00062 { 00063 ErrorToPut = 0.025f; 00064 } 00065 else 00066 { 00067 ErrorToPut = 0.05f; 00068 } 00069 iso_put_error(ipISO,nelem,iso.numLevels_max[ipISO][nelem],ipHi,IPRAD,ErrorToPut); 00070 } 00071 /* Error for total recombination */ 00072 iso_put_error(ipISO,nelem,iso.numLevels_max[ipISO][nelem],iso.numLevels_max[ipISO][nelem],IPRAD,0.02f); 00073 00074 for( ipHi=1; ipHi<= iso.numLevels_max[ipISO][nelem]; ipHi++ ) 00075 { 00076 /* >>chng 06 mar 15, the upper limit incorrectly went to numLevels_max */ 00077 for( ipLo=0; ipLo<= ipHi; ipLo++ ) 00078 { 00079 for( typeOfRate=0; typeOfRate<=1; typeOfRate++ ) 00080 { 00081 if( iso.Error[ipISO][nelem][ipHi][ipLo][typeOfRate] >= 0. ) 00082 { 00083 iso.ErrorFactor[ipISO][nelem][ipHi][ipLo][typeOfRate] = 00084 (realnum)MyGaussRand( iso.Error[ipISO][nelem][ipHi][ipLo][typeOfRate] ); 00085 ASSERT( iso.ErrorFactor[ipISO][nelem][ipHi][ipLo][typeOfRate] > 0. ); 00086 } 00087 else 00088 { 00089 iso.ErrorFactor[ipISO][nelem][ipHi][ipLo][typeOfRate] = 1.0f; 00090 } 00091 } 00092 } 00093 } 00094 00095 /* set flag saying that error generation has been done. */ 00096 iso.lgErrGenDone[ipISO][nelem] = true; 00097 return; 00098 }