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 /*ipoint returns pointer to any energy within energy mesh */ 00004 /*ipFineCont returns array index within fine energy mesh */ 00005 /*ipContEnergy generate unique pointer to energy within continuum array 00006 * continuum energy in Rydbergs */ 00007 /*ipLineEnergy generate unique pointer to line energy within energy mesh 00008 * line energy in Rydbergs */ 00009 #include "cddefines.h" 00010 #include "continuum.h" 00011 #include "prt.h" 00012 #include "rfield.h" 00013 #include "ipoint.h" 00014 00015 /*ipoint returns pointer to any energy within energy mesh */ 00016 long ipoint(double energy_ryd) 00017 { 00018 long int i, 00019 ipoint_v; 00020 00021 DEBUG_ENTRY( "ipoint()" ); 00022 00023 if( energy_ryd < continuum.filbnd[0] || energy_ryd > continuum.filbnd[continuum.nrange] ) 00024 { 00025 fprintf( ioQQQ, " ipoint:\n" ); 00026 fprintf( ioQQQ, " The energy_ryd array is not defined at nu=%11.3e. The bounds are%11.3e%11.3e\n", 00027 energy_ryd, continuum.filbnd[0], continuum.filbnd[continuum.nrange] ); 00028 fprintf( ioQQQ, " ipoint is aborting to get trace, to find how this happened\n" ); 00029 ShowMe(); 00030 cdEXIT(EXIT_FAILURE); 00031 } 00032 00033 for( i=0; i < continuum.nrange; i++ ) 00034 { 00035 if( energy_ryd >= continuum.filbnd[i] && energy_ryd <= continuum.filbnd[i+1] ) 00036 { 00037 00038 /* this is on the Fortran scale of array index counting, so is one above the 00039 * c scale. later the -1 will appear on any index references */ 00040 ipoint_v = (long int)(log10(energy_ryd/continuum.filbnd[i])/continuum.fildel[i] + 00041 1.0 + continuum.ifill0[i]); 00042 00043 ASSERT( ipoint_v >= 0 ); 00044 /* recall on F scale */ 00045 ipoint_v = MIN2( rfield.nupper , ipoint_v ); 00046 return ipoint_v; 00047 } 00048 } 00049 00050 /* this exit not possible, here to shut up some compilers */ 00051 fprintf( ioQQQ, " IPOINT logic error, energy=%.2e\n", 00052 energy_ryd ); 00053 cdEXIT(EXIT_FAILURE); 00054 } 00055 00056 /*ipContEnergy generate unique pointer to energy within continuum array */ 00057 long ipContEnergy( 00058 /* continuum energy in Rydbergs */ 00059 double energy, 00060 /* 4 char label for continuum, like those returned by chLineLbl */ 00061 const char *chLabel) 00062 { 00063 long int ipConSafe_v; 00064 00065 DEBUG_ENTRY( "ipContEnergy()" ); 00066 00067 ipConSafe_v = ipoint(energy); 00068 00069 /* write label in this cell if not anything there yet */ 00070 if( strcmp(rfield.chContLabel[ipConSafe_v-1]," ") == 0 ) 00071 { 00072 strcpy( rfield.chContLabel[ipConSafe_v-1], chLabel ); 00073 } 00074 00075 /* this is a quick way to find all continua that occur within a given cell, 00076 * recall the off-by-one aspect of C */ 00077 { 00078 enum {DEBUG_LOC=false}; 00079 if( DEBUG_LOC ) 00080 { 00081 /* recall the off-by-one aspect of C - the number below is 00082 * on the physics scale because this returns the index 00083 * on that scale, so the first is 1 for [0] */ 00084 if( ipConSafe_v == 23 ) 00085 fprintf(ioQQQ,"%s\n", chLabel ); 00086 } 00087 } 00088 return ipConSafe_v; 00089 } 00090 00091 /*ipLineEnergy generate unique pointer to line energy within energy mesh 00092 * line energy in Rydbergs - 00093 * return value is array index on the physics or Fortran scale, counting from 1 */ 00094 long ipLineEnergy(double energy, 00095 /* 4 char label for line, like those returned by chLineLbl */ 00096 const char *chLabel , 00097 /* will make sure energy is < this array index if greater than 0, does nothing if <= 0*/ 00098 long ipIonEnergy ) 00099 { 00100 long int ipLine_ret; 00101 00102 DEBUG_ENTRY( "ipLineEnergy()" ); 00103 00104 ipLine_ret = ipoint(energy); 00105 ASSERT( ipLine_ret ); 00106 /* make sure pointer is below next higher continuum if positive */ 00107 if( ipIonEnergy > 0 ) 00108 { 00109 ipLine_ret = MIN2( ipLine_ret , ipIonEnergy-1 ); 00110 } 00111 00112 ASSERT( ipLine_ret > 0 ); 00113 /* stuff in a label if none there, 00114 * note that this is offset one below actual number to be on C scale of arrays */ 00115 /* >>chng 06 nov 23, faster to use line_count index rather than checking 5 chars 00116 * first call will have zero so false */ 00117 /*if( strcmp(rfield.chLineLabel[ipLine_ret-1]," ")==0 )*/ 00118 if( !rfield.line_count[ipLine_ret-1] ) 00119 { 00120 strcpy( rfield.chLineLabel[ipLine_ret-1], chLabel ); 00121 } 00122 /* this keeps track of the number of lines within this cell */ 00123 ++rfield.line_count[ipLine_ret-1]; 00124 00125 /* this is a quick way to find all lines that occur within a given cell, 00126 * recall the off-by-one aspect of C */ 00127 { 00128 enum {DEBUG_LOC=false}; 00129 if( DEBUG_LOC ) 00130 { 00131 /* recall the off-by-one aspect of C - the numbers is one more 00132 * than the index on the C scale */ 00133 if( ipLine_ret == 23 ) 00134 fprintf(ioQQQ,"%s\n", chLabel ); 00135 } 00136 } 00137 00138 /* this implements the print continuum indices command */ 00139 if( prt.lgPrtContIndices ) 00140 { 00141 /* print header if first time */ 00142 static bool lgFirst = true; 00143 if( lgFirst ) 00144 { 00145 /* print header and set flag so do not do again */ 00146 fprintf(ioQQQ , "\n\noutput from print continuum indices command follows.\n"); 00147 fprintf(ioQQQ , "cont ind (F scale)\tenergy(ryd)\tlabel\n"); 00148 lgFirst = false; 00149 } 00150 if( energy >= prt.lgPrtContIndices_lo_E && energy <= prt.lgPrtContIndices_hi_E ) 00151 { 00152 /* use varying formats depending on order of magnitude of energy 00153 * NB - the ipLine_ret is the index on the physical or Fortran scale, 00154 * and is 1 greater than the array element used in the code, which is 00155 * on the C scale */ 00156 if( energy < 1. ) 00157 { 00158 fprintf(ioQQQ , "%li\t%.3e\t%s\n" , ipLine_ret , energy , chLabel); 00159 } 00160 else if( energy < 10. ) 00161 { 00162 fprintf(ioQQQ , "%li\t%.3f\t%s\n" , ipLine_ret , energy , chLabel); 00163 } 00164 else if( energy < 100. ) 00165 { 00166 fprintf(ioQQQ , "%li\t%.2f\t%s\n" , ipLine_ret , energy , chLabel); 00167 } 00168 else 00169 { 00170 fprintf(ioQQQ , "%li\t%.1f\t%s\n" , ipLine_ret , energy , chLabel); 00171 } 00172 } 00173 } 00174 00175 if( prt.lgPrnLineCell ) 00176 { 00177 /* print line cell option - number on command line is cell on Physics scale */ 00178 if( prt.nPrnLineCell == ipLine_ret ) 00179 { 00180 static bool lgMustPrintHeader = true; 00181 if( lgMustPrintHeader ) 00182 fprintf(ioQQQ, "Lines within cell %li\n",prt.nPrnLineCell ); 00183 lgMustPrintHeader = false; 00184 fprintf(ioQQQ,"**>%s<**\n" , chLabel ); 00185 } 00186 } 00187 return ipLine_ret; 00188 } 00189 00190 /*ipFineCont returns array index within fine energy mesh */ 00191 long ipFineCont( 00192 /* energy in Ryd */ 00193 double energy_ryd ) 00194 { 00195 long int ipoint_v; 00196 00197 DEBUG_ENTRY( "ipFineCont()" ); 00198 00199 if( energy_ryd < rfield.fine_ener_lo || energy_ryd > rfield.fine_ener_hi ) 00200 { 00201 return -1; 00202 } 00203 00204 /* this is on the Fortran scale of array index counting, so is one above the 00205 * c scale. later the -1 will appear on any index references 00206 * 00207 * 07 Jun 22 testing done to confirm that energy grid is correct: did 00208 * same sim with standard fine continuum resolution, and another with 200x 00209 * finer resolution, and confirmed that these line up correctly. */ 00210 ipoint_v = (long int)(log10(energy_ryd*(1.-rfield.fine_resol/2.) / 00211 rfield.fine_ener_lo)/log10(1.+rfield.fine_resol)); 00212 00213 ASSERT( ipoint_v >= 0 && ipoint_v< rfield.nfine_malloc ); 00214 return ipoint_v; 00215 }