cloudy trunk
Loading...
Searching...
No Matches
cont_ipoint.cpp
Go to the documentation of this file.
1/* This file is part of Cloudy and is copyright (C)1978-2013 by Gary J. Ferland and
2 * others. For conditions of distribution and use see copyright notice in license.txt */
3/*ipoint returns pointer to any energy within energy mesh */
4/*ipFine()Cont returns array index within fine energy mesh */
5/*ipContEnergy generate pointer to energy within continuum array
6 * continuum energy in Rydbergs */
7/*ipLineEnergy generate pointer to line energy within energy mesh
8 * line energy in Rydbergs */
9#include "cddefines.h"
10#include "continuum.h"
11#include "prt.h"
12#include "rfield.h"
13#include "ipoint.h"
14
15/*ipoint returns pointer to any energy within energy mesh */
16long ipoint(double energy_ryd)
17{
18 long int i,
19 ipoint_v;
20
21 DEBUG_ENTRY( "ipoint()" );
22
23 // make sure mesh is set up
24 ASSERT( continuum.nrange > 0 );
25
26 if( energy_ryd < continuum.filbnd[0] || energy_ryd > continuum.filbnd[continuum.nrange] )
27 {
28 fprintf( ioQQQ, " ipoint:\n" );
29 fprintf( ioQQQ, " The energy_ryd array is not defined at nu=%11.3e. The bounds are%11.3e%11.3e\n",
30 energy_ryd, continuum.filbnd[0], continuum.filbnd[continuum.nrange] );
31 fprintf( ioQQQ, " ipoint is aborting to get trace, to find how this happened\n" );
32 ShowMe();
34 }
35
36 for( i=0; i < continuum.nrange; i++ )
37 {
38 if( energy_ryd >= continuum.filbnd[i] && energy_ryd <= continuum.filbnd[i+1] )
39 {
40
41 /* this is on the Fortran scale of array index counting, so is one above the
42 * c scale. later the -1 will appear on any index references */
43 ipoint_v = (long int)(log10(energy_ryd/continuum.filbnd[i])/continuum.fildel[i] +
44 1.0 + continuum.ifill0[i]);
45
46 ASSERT( ipoint_v >= 0 );
47 /* recall on F scale */
48 ipoint_v = MIN2( rfield.nupper , ipoint_v );
49
50 if( ipoint_v < rfield.nflux-2 && ipoint_v>2 )
51 {
52 // possible will need to adjust index if cells have been fiddled with
53 if( energy_ryd > rfield.anu[ipoint_v-1]+rfield.widflx[ipoint_v-1]/2. )
54 ++ipoint_v;
55 if( energy_ryd < rfield.anu[ipoint_v-1]-rfield.widflx[ipoint_v-1]/2. )
56 --ipoint_v;
57 {
58 enum {DEBUG_LOC=false};
59 if( DEBUG_LOC )
60 {
61 if( energy_ryd > rfield.anu[ipoint_v]+rfield.widflx[ipoint_v]/2. )
62 {
63 fprintf(ioQQQ,"DISASTER ipoint hi bounds about to throw "\
64 "energy_ryd=%e hi bound=%e ipoint_v=%li nflux=%li\n",
65 energy_ryd , rfield.anu[ipoint_v]+rfield.widflx[ipoint_v]/2.,
66 ipoint_v,rfield.nflux);
67 }
68 if( energy_ryd < rfield.anu[ipoint_v-2]-rfield.widflx[ipoint_v-2]/2. )
69 {
70 fprintf(ioQQQ,"DISASTER ipoint low bounds about to throw "\
71 "energy_ryd=%e low bound=%e ipoint_v=%li nflux=%li\n",
72 energy_ryd , rfield.anu[ipoint_v-2]-rfield.widflx[ipoint_v-2]/2.,
73 ipoint_v,rfield.nflux);
74 }
75 }
76 }
77
78 ASSERT( energy_ryd <= rfield.anu[ipoint_v]+rfield.widflx[ipoint_v]/2. );
79 ASSERT( energy_ryd >= rfield.anu[ipoint_v-2]-rfield.widflx[ipoint_v-2]/2. );
80 }
81 return ipoint_v;
82 }
83 }
84
85 /* this exit not possible, here to shut up some compilers */
86 fprintf( ioQQQ, " IPOINT logic error, energy=%.2e\n",
87 energy_ryd );
89}
90
91/*ipContEnergy generate pointer to energy within continuum array */
93 /* continuum energy in Rydbergs */
94 double energy,
95 /* 4 char label for continuum, like those returned by chLineLbl */
96 const char *chLabel)
97{
98 long int ipConSafe_v;
99
100 DEBUG_ENTRY( "ipContEnergy()" );
101
102 ipConSafe_v = ipoint(energy);
103
104 /* write label in this cell if not anything there yet */
105 if( strcmp(rfield.chContLabel[ipConSafe_v-1]," ") == 0 )
106 {
107 strcpy( rfield.chContLabel[ipConSafe_v-1], chLabel );
108 }
109
110 /* this is a quick way to find all continua that occur within a given cell,
111 * recall the off-by-one aspect of C */
112 {
113 enum {DEBUG_LOC=false};
114 if( DEBUG_LOC )
115 {
116 /* recall the off-by-one aspect of C - the number below is
117 * on the physics scale because this returns the index
118 * on that scale, so the first is 1 for [0] */
119 if( ipConSafe_v == 23 )
120 fprintf(ioQQQ,"%s\n", chLabel );
121 }
122 }
123 return ipConSafe_v;
124}
125
126/*ipLineEnergy generate pointer to line energy within energy mesh
127 * line energy in Rydbergs -
128 * return value is array index on the physics or Fortran scale, counting from 1 */
129long ipLineEnergy(double energy,
130 /* 4 char label for line, like those returned by chLineLbl */
131 const char *chLabel ,
132 /* will make sure energy is < this array index if greater than 0, does nothing if <= 0*/
133 long ipIonEnergy )
134{
135 long int ipLine_ret;
136
137 DEBUG_ENTRY( "ipLineEnergy()" );
138
139 ipLine_ret = ipoint(energy);
140 ASSERT( ipLine_ret );
141 /* make sure pointer is below next higher continuum if positive */
142 if( ipIonEnergy > 0 )
143 {
144 ipLine_ret = MIN2( ipLine_ret , ipIonEnergy-1 );
145 }
146
147 ASSERT( ipLine_ret > 0 );
148 /* stuff in a label if none there,
149 * note that this is offset one below actual number to be on C scale of arrays */
150 /* >>chng 06 nov 23, faster to use line_count index rather than checking 5 chars
151 * first call will have zero so false */
152 /*if( strcmp(rfield.chLineLabel[ipLine_ret-1]," ")==0 )*/
153 if( !rfield.line_count[ipLine_ret-1] )
154 {
155 strcpy( rfield.chLineLabel[ipLine_ret-1], chLabel );
156 }
157 /* this keeps track of the number of lines within this cell */
158 ++rfield.line_count[ipLine_ret-1];
159
160 /* this is a quick way to find all lines that occur within a given cell,
161 * recall the off-by-one aspect of C */
162 {
163 enum {DEBUG_LOC=false};
164 if( DEBUG_LOC )
165 {
166 /* recall the off-by-one aspect of C - the numbers is one more
167 * than the index on the C scale */
168 if( ipLine_ret == 23 )
169 fprintf(ioQQQ,"%s\n", chLabel );
170 }
171 }
172
173 /* this implements the print continuum indices command */
174 if( prt.lgPrtContIndices )
175 {
176 /* print header if first time */
177 static bool lgFirst = true;
178 if( lgFirst )
179 {
180 /* print header and set flag so do not do again */
181 fprintf(ioQQQ , "\n\noutput from print continuum indices command follows.\n");
182 fprintf(ioQQQ , "cont ind (F scale)\tenergy(ryd)\tlabel\n");
183 lgFirst = false;
184 }
185 if( energy >= prt.lgPrtContIndices_lo_E && energy <= prt.lgPrtContIndices_hi_E )
186 {
187 /* use varying formats depending on order of magnitude of energy
188 * NB - the ipLine_ret is the index on the physical or Fortran scale,
189 * and is 1 greater than the array element used in the code, which is
190 * on the C scale */
191 if( energy < 1. )
192 {
193 fprintf(ioQQQ , "%li\t%.3e\t%s\n" , ipLine_ret , energy , chLabel);
194 }
195 else if( energy < 10. )
196 {
197 fprintf(ioQQQ , "%li\t%.3f\t%s\n" , ipLine_ret , energy , chLabel);
198 }
199 else if( energy < 100. )
200 {
201 fprintf(ioQQQ , "%li\t%.2f\t%s\n" , ipLine_ret , energy , chLabel);
202 }
203 else
204 {
205 fprintf(ioQQQ , "%li\t%.1f\t%s\n" , ipLine_ret , energy , chLabel);
206 }
207 }
208 }
209
210 if( prt.lgPrnLineCell )
211 {
212 /* print line cell option - number on command line is cell on Physics scale */
213 if( prt.nPrnLineCell == ipLine_ret )
214 {
215 static bool lgMustPrintHeader = true;
217 fprintf(ioQQQ, "Lines within cell %li (physics scale) \nLabel\tEnergy(Ryd)\n",prt.nPrnLineCell );
218 lgMustPrintHeader = false;
219 fprintf(ioQQQ,"%s\t%.3e\n" , chLabel , energy );
220 }
221 }
222 return ipLine_ret;
223}
224
225/*ipFine()Cont returns array index within fine energy mesh */
227 /* energy in Ryd */
228 double energy_ryd )
229{
230 long int ipoint_v;
231
232 DEBUG_ENTRY( "ipFine()Cont()" );
233
234 if( energy_ryd < rfield.fine_ener_lo || energy_ryd > rfield.fine_ener_hi )
235 {
236 return -1;
237 }
238
239 /* this is on the Fortran scale of array index counting, so is one above the
240 * c scale. later the -1 will appear on any index references
241 *
242 * 07 Jun 22 testing done to confirm that energy grid is correct: did
243 * same sim with standard fine continuum resolution, and another with 200x
244 * finer resolution, and confirmed that these line up correctly. */
245 ipoint_v = (long int)(log10(energy_ryd*(1.-rfield.fine_resol/2.) /
246 rfield.fine_ener_lo)/log10(1.+rfield.fine_resol));
247
248 ASSERT( ipoint_v >= 0 && ipoint_v< rfield.nfine_malloc );
249 return ipoint_v;
250}
FILE * ioQQQ
Definition cddefines.cpp:7
#define ASSERT(exp)
Definition cddefines.h:578
#define MIN2
Definition cddefines.h:761
#define EXIT_FAILURE
Definition cddefines.h:140
#define cdEXIT(FAIL)
Definition cddefines.h:434
#define DEBUG_ENTRY(funcname)
Definition cddefines.h:684
void ShowMe(void)
Definition service.cpp:181
long ipFineCont(double energy_ryd)
long ipoint(double energy_ryd)
long ipLineEnergy(double energy, const char *chLabel, long ipIonEnergy)
long ipContEnergy(double energy, const char *chLabel)
t_continuum continuum
Definition continuum.cpp:5
t_prt prt
Definition prt.cpp:10
static bool lgFirst
t_rfield rfield
Definition rfield.cpp:8
static bool lgMustPrintHeader