00001
00002
00003
00004
00005
00006
00007
00008
00009 #include "cddefines.h"
00010 #include "continuum.h"
00011 #include "prt.h"
00012 #include "rfield.h"
00013 #include "ipoint.h"
00014
00015
00016 long ipoint(double energy_ryd)
00017 {
00018 long int i,
00019 ipoint_v;
00020
00021 DEBUG_ENTRY( "ipoint()" );
00022
00023
00024 ASSERT( continuum.nrange > 0 );
00025
00026 if( energy_ryd < continuum.filbnd[0] || energy_ryd > continuum.filbnd[continuum.nrange] )
00027 {
00028 fprintf( ioQQQ, " ipoint:\n" );
00029 fprintf( ioQQQ, " The energy_ryd array is not defined at nu=%11.3e. The bounds are%11.3e%11.3e\n",
00030 energy_ryd, continuum.filbnd[0], continuum.filbnd[continuum.nrange] );
00031 fprintf( ioQQQ, " ipoint is aborting to get trace, to find how this happened\n" );
00032 ShowMe();
00033 cdEXIT(EXIT_FAILURE);
00034 }
00035
00036 for( i=0; i < continuum.nrange; i++ )
00037 {
00038 if( energy_ryd >= continuum.filbnd[i] && energy_ryd <= continuum.filbnd[i+1] )
00039 {
00040
00041
00042
00043 ipoint_v = (long int)(log10(energy_ryd/continuum.filbnd[i])/continuum.fildel[i] +
00044 1.0 + continuum.ifill0[i]);
00045
00046 ASSERT( ipoint_v >= 0 );
00047
00048 ipoint_v = MIN2( rfield.nupper , ipoint_v );
00049
00050 if( ipoint_v < rfield.nflux-2 && ipoint_v>2 )
00051 {
00052
00053 if( energy_ryd > rfield.anu[ipoint_v-1]+rfield.widflx[ipoint_v-1]/2. )
00054 ++ipoint_v;
00055 if( energy_ryd < rfield.anu[ipoint_v-1]-rfield.widflx[ipoint_v-1]/2. )
00056 --ipoint_v;
00057 {
00058 enum {DEBUG_LOC=false};
00059 if( DEBUG_LOC )
00060 {
00061 if( energy_ryd > rfield.anu[ipoint_v]+rfield.widflx[ipoint_v]/2. )
00062 {
00063 fprintf(ioQQQ,"DISASTER ipoint hi bounds about to throw "\
00064 "energy_ryd=%e hi bound=%e ipoint_v=%li nflux=%li\n",
00065 energy_ryd , rfield.anu[ipoint_v]+rfield.widflx[ipoint_v]/2.,
00066 ipoint_v,rfield.nflux);
00067 }
00068 if( energy_ryd < rfield.anu[ipoint_v-2]-rfield.widflx[ipoint_v-2]/2. )
00069 {
00070 fprintf(ioQQQ,"DISASTER ipoint low bounds about to throw "\
00071 "energy_ryd=%e low bound=%e ipoint_v=%li nflux=%li\n",
00072 energy_ryd , rfield.anu[ipoint_v-2]-rfield.widflx[ipoint_v-2]/2.,
00073 ipoint_v,rfield.nflux);
00074 }
00075 }
00076 }
00077
00078 ASSERT( energy_ryd <= rfield.anu[ipoint_v]+rfield.widflx[ipoint_v]/2. );
00079 ASSERT( energy_ryd >= rfield.anu[ipoint_v-2]-rfield.widflx[ipoint_v-2]/2. );
00080 }
00081 return ipoint_v;
00082 }
00083 }
00084
00085
00086 fprintf( ioQQQ, " IPOINT logic error, energy=%.2e\n",
00087 energy_ryd );
00088 cdEXIT(EXIT_FAILURE);
00089 }
00090
00091
00092 long ipContEnergy(
00093
00094 double energy,
00095
00096 const char *chLabel)
00097 {
00098 long int ipConSafe_v;
00099
00100 DEBUG_ENTRY( "ipContEnergy()" );
00101
00102 ipConSafe_v = ipoint(energy);
00103
00104
00105 if( strcmp(rfield.chContLabel[ipConSafe_v-1]," ") == 0 )
00106 {
00107 strcpy( rfield.chContLabel[ipConSafe_v-1], chLabel );
00108 }
00109
00110
00111
00112 {
00113 enum {DEBUG_LOC=false};
00114 if( DEBUG_LOC )
00115 {
00116
00117
00118
00119 if( ipConSafe_v == 23 )
00120 fprintf(ioQQQ,"%s\n", chLabel );
00121 }
00122 }
00123 return ipConSafe_v;
00124 }
00125
00126
00127
00128
00129 long ipLineEnergy(double energy,
00130
00131 const char *chLabel ,
00132
00133 long ipIonEnergy )
00134 {
00135 long int ipLine_ret;
00136
00137 DEBUG_ENTRY( "ipLineEnergy()" );
00138
00139 ipLine_ret = ipoint(energy);
00140 ASSERT( ipLine_ret );
00141
00142 if( ipIonEnergy > 0 )
00143 {
00144 ipLine_ret = MIN2( ipLine_ret , ipIonEnergy-1 );
00145 }
00146
00147 ASSERT( ipLine_ret > 0 );
00148
00149
00150
00151
00152
00153 if( !rfield.line_count[ipLine_ret-1] )
00154 {
00155 strcpy( rfield.chLineLabel[ipLine_ret-1], chLabel );
00156 }
00157
00158 ++rfield.line_count[ipLine_ret-1];
00159
00160
00161
00162 {
00163 enum {DEBUG_LOC=false};
00164 if( DEBUG_LOC )
00165 {
00166
00167
00168 if( ipLine_ret == 23 )
00169 fprintf(ioQQQ,"%s\n", chLabel );
00170 }
00171 }
00172
00173
00174 if( prt.lgPrtContIndices )
00175 {
00176
00177 static bool lgFirst = true;
00178 if( lgFirst )
00179 {
00180
00181 fprintf(ioQQQ , "\n\noutput from print continuum indices command follows.\n");
00182 fprintf(ioQQQ , "cont ind (F scale)\tenergy(ryd)\tlabel\n");
00183 lgFirst = false;
00184 }
00185 if( energy >= prt.lgPrtContIndices_lo_E && energy <= prt.lgPrtContIndices_hi_E )
00186 {
00187
00188
00189
00190
00191 if( energy < 1. )
00192 {
00193 fprintf(ioQQQ , "%li\t%.3e\t%s\n" , ipLine_ret , energy , chLabel);
00194 }
00195 else if( energy < 10. )
00196 {
00197 fprintf(ioQQQ , "%li\t%.3f\t%s\n" , ipLine_ret , energy , chLabel);
00198 }
00199 else if( energy < 100. )
00200 {
00201 fprintf(ioQQQ , "%li\t%.2f\t%s\n" , ipLine_ret , energy , chLabel);
00202 }
00203 else
00204 {
00205 fprintf(ioQQQ , "%li\t%.1f\t%s\n" , ipLine_ret , energy , chLabel);
00206 }
00207 }
00208 }
00209
00210 if( prt.lgPrnLineCell )
00211 {
00212
00213 if( prt.nPrnLineCell == ipLine_ret )
00214 {
00215 static bool lgMustPrintHeader = true;
00216 if( lgMustPrintHeader )
00217 fprintf(ioQQQ, "Lines within cell %li (physics scale) \nLabel\tEnergy(Ryd)\n",prt.nPrnLineCell );
00218 lgMustPrintHeader = false;
00219 fprintf(ioQQQ,"%s\t%.3e\n" , chLabel , energy );
00220 }
00221 }
00222 return ipLine_ret;
00223 }
00224
00225
00226 long ipFineCont(
00227
00228 double energy_ryd )
00229 {
00230 long int ipoint_v;
00231
00232 DEBUG_ENTRY( "ipFine()Cont()" );
00233
00234 if( energy_ryd < rfield.fine_ener_lo || energy_ryd > rfield.fine_ener_hi )
00235 {
00236 return -1;
00237 }
00238
00239
00240
00241
00242
00243
00244
00245 ipoint_v = (long int)(log10(energy_ryd*(1.-rfield.fine_resol/2.) /
00246 rfield.fine_ener_lo)/log10(1.+rfield.fine_resol));
00247
00248 ASSERT( ipoint_v >= 0 && ipoint_v< rfield.nfine_malloc );
00249 return ipoint_v;
00250 }