physics.cpp 12.3 KB
Newer Older
1

steffen.schotthoefer's avatar
steffen.schotthoefer committed
2
#include "physics.h"
steffen.schotthoefer's avatar
steffen.schotthoefer committed
3
#include "spline.h"
4
#include "toolboxes/errormessages.h"
steffen.schotthoefer's avatar
steffen.schotthoefer committed
5

steffen.schotthoefer's avatar
steffen.schotthoefer committed
6
#include <fstream>
steffen.schotthoefer's avatar
steffen.schotthoefer committed
7
#include <list>
steffen.schotthoefer's avatar
steffen.schotthoefer committed
8
#include <map>
9

10
11
Physics::Physics() {
    // LoadDatabase
12
13
}

14
void Physics::LoadDatabase( std::string fileName_H, std::string fileName_O, std::string fileName_stppower ) {
15
16
    VectorVector transport_XS_H;
    VectorVector transport_XS_O;
steffen.schotthoefer's avatar
steffen.schotthoefer committed
17

18
    VectorVector scattering_XS_H;
19
    VectorVector scattering_XS_O;
20
21

    VectorVector total_scat_XS_H;
22
    VectorVector total_scat_XS_O;
23

jannick.wolters's avatar
jannick.wolters committed
24
25
26
27
    std::vector<VectorVector> headers_H;
    std::vector<VectorVector> data_H;
    std::vector<VectorVector> headers_O;
    std::vector<VectorVector> data_O;
28
29

    // Read data for H
jannick.wolters's avatar
jannick.wolters committed
30
    std::tie( headers_H, data_H ) = ReadENDL( fileName_H );
31

jannick.wolters's avatar
jannick.wolters committed
32
33
34
    for( unsigned i = 0; i < headers_H.size(); ++i ) {
        auto header = headers_H[i];
        auto data   = data_H[i];
35
36
37
38
39
40
41
42
43
44
45
46
47
48
        if( header[1][0] == 7 && header[1][1] == 0 ) {
            transport_XS_H = data;
        }
        // Integrated elastic scattering XS
        else if( header[1][0] == 10 && header[1][1] == 0 ) {
            total_scat_XS_H = data;
        }
        // Angular distribution large angle elastic scattering XS
        else if( header[1][0] == 8 && header[1][1] == 22 ) {
            scattering_XS_H = data;
        }
    }

    // Read data for O
jannick.wolters's avatar
jannick.wolters committed
49
    std::tie( headers_O, data_O ) = ReadENDL( fileName_O );
50
51

    // Find required quantities and transfer to matrix
52
    for( unsigned i = 0; i < headers_O.size(); ++i ) {
jannick.wolters's avatar
jannick.wolters committed
53
54
        auto header = headers_O[i];
        auto data   = data_O[i];
55
        if( header[1][0] == 7 && header[1][1] == 0 ) {
56
            transport_XS_O = data;
57
58
59
        }
        // Integrated elastic scattering XS
        else if( header[1][0] == 10 && header[1][1] == 0 ) {
60
            total_scat_XS_O = data;
61
62
63
        }
        // Angular distribution large angle elastic scattering XS
        else if( header[1][0] == 8 && header[1][1] == 22 ) {
64
            scattering_XS_O = data;
65
66
67
68
69
70
71
72
73
        }
    }

    if( transport_XS_H.size() != transport_XS_O.size() )
        ErrorMessages::Error( "transport_XS of hydrogen and oxygen have different sizes", CURRENT_FUNCTION );
    if( scattering_XS_H.size() != scattering_XS_O.size() )
        ErrorMessages::Error( "scattering_XS of hydrogen and oxygen have different sizes", CURRENT_FUNCTION );
    if( total_scat_XS_H.size() != total_scat_XS_O.size() )
        ErrorMessages::Error( "total_XS of hydrogen and oxygen have different sizes", CURRENT_FUNCTION );
74

75
76
77
78
79
    // Combine values for H and O according to mixture ratio in water
    for( unsigned i = 0; i < transport_XS_H.size(); ++i )
        _xsTransportH2O.push_back( 0.11189400 * transport_XS_H[i] + 0.88810600 * transport_XS_O[i] );
    for( unsigned i = 0; i < scattering_XS_H.size(); ++i ) _xsH2O.push_back( 0.11189400 * scattering_XS_H[i] + 0.88810600 * scattering_XS_O[i] );
    for( unsigned i = 0; i < total_scat_XS_H.size(); ++i ) _xsTotalH2O.push_back( 0.11189400 * total_scat_XS_H[i] + 0.88810600 * total_scat_XS_O[i] );
80
    std::cout << "size is " << _xsH2O.size() << " " << _xsH2O[0].size() << std::endl;
81
    _stpowH2O = ReadStoppingPowers( fileName_stppower );
82
83
}

84
VectorVector Physics::GetScatteringXS( Vector energies, Vector angle ) {
85

86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
    std::cout << _xsH2O.size() << " " << _xsH2O[0].size() << std::endl;
    std::vector<std::vector<double>> tmp;              // vectorvector which stores data at fixed energies
    std::vector<std::vector<double>> xsH2OGrid;        // matrix which stores tensorized data for given angular grid, original energy grid
    std::vector<std::vector<double>> xsH2OGridGrid;    // matrix which stores tensorized data for given angular and energy grid
    std::vector<double> tmpAngleGrid;
    std::vector<double> tmpEnergyGrid;
    std::vector<double> tmp1;
    std::vector<double> energiesOrig;    // original energy grid
    tmpAngleGrid.resize( angle.size() );
    tmp.clear();
    tmp.resize( 0 );
    xsH2OGrid.resize( 0 );
    double energyCurrent = _xsH2O[0][0];

    // build grid at original energies and given angular grid
    for( unsigned i = 0; i < _xsH2O.size(); ++i ) {
        // split vector into subvectors with identical energy
        if( abs( _xsH2O[i][0] - energyCurrent ) < 1e-12 ) {
            tmp.push_back( std::vector<double>( { _xsH2O[i][1], _xsH2O[i][2] } ) );
        }
        else {
            // new energy section starts in _xsH2O. Now we interpolate the values in tmp at given angular grid

            // interpolate vector at different angles for fixed current energies
            Spline interp;
            interp.set_points( tmp[0], tmp[1], false );    // false == linear interpolation
            for( unsigned k = 0; k < angle.size(); k++ ) {
                // why do we need if else statements?
                // Linear interpolation
                tmpAngleGrid[k] = interp( energies[k] );
            }
            xsH2OGrid.push_back( tmpAngleGrid );

            // reset current energy
            energiesOrig.push_back( energyCurrent );
            energyCurrent = _xsH2O[i][0];
            tmp.clear();
            tmp.resize( 0 );
            tmp.push_back( std::vector<double>( { _xsH2O[i][1], _xsH2O[i][2] } ) );
        }
    }
    // perform interpolation at fixed original energy for new energy grid

    for( unsigned j = 0; j < angle.size(); ++j ) {    // loop over all angles
        tmp1.resize( 0 );
        // store all values at given angle j for all original energies
        for( unsigned i = 0; i < energiesOrig.size(); ++i ) tmp1.push_back( xsH2OGrid[i][j] );
        Spline interp;
        interp.set_points( energiesOrig, tmp1, false );    // false == linear interpolation
        for( unsigned k = 0; k < energies.size(); k++ ) {
            // why do we need if else statements?
            // Linear interpolation
            tmpEnergyGrid[k] = interp( energies[k] );
        }
        xsH2OGridGrid.push_back( tmpEnergyGrid );
    }

    //_xsH2O
144
145
146
147
148
149
150
151
152
    VectorVector scattering_XS( xsH2OGridGrid.size() );
    // write vector<vector> to VectorVector
    for( unsigned idx_energy = 0; idx_energy < xsH2OGridGrid.size(); idx_energy++ ) {
        scattering_XS[idx_energy] = Vector( xsH2OGridGrid[idx_energy].size() );
        for( unsigned idx_angular = 0; idx_angular < xsH2OGridGrid[idx_energy].size(); idx_angular++ ) {
            scattering_XS[idx_energy][idx_angular] = xsH2OGridGrid[idx_energy][idx_angular];
        }
    }

pia.stammer's avatar
pia.stammer committed
153
    return scattering_XS;
154
155
}

156
VectorVector Physics::GetTotalXS( Vector energies, Vector density ) {
pia.stammer's avatar
pia.stammer committed
157
    VectorVector total_XS;
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
    double xsH2O_i;

    Spline interp;
    interp.set_points( _xsTotalH2O[0], _xsTotalH2O[1], false );    // false == linear interpolation

    for( unsigned i = 0; i < energies.size(); i++ ) {

        if( energies[i] < _xsTotalH2O[1][0] ) {
            xsH2O_i = _xsTotalH2O[1][0];
        }
        else if( energies[i] > _xsTotalH2O[1][energies.size() - 1] ) {
            xsH2O_i = _xsTotalH2O[1][energies.size() - 1];
        }
        else {
            // Linear interpolation
            xsH2O_i = interp( energies[i] );
        }

        total_XS[i] = xsH2O_i * density;
    }

pia.stammer's avatar
pia.stammer committed
179
    return total_XS;
180
181
}

182
VectorVector Physics::GetStoppingPower( Vector energies, Vector density ) {
pia.stammer's avatar
pia.stammer committed
183
    VectorVector stopping_power;
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
    double stpw_H2O_i;

    Spline interp;
    interp.set_points( _stpowH2O[0], _stpowH2O[1], false );    // false == linear interpolation

    for( unsigned i = 0; i < energies.size(); i++ ) {

        if( energies[i] < _stpowH2O[1][0] ) {
            stpw_H2O_i = _stpowH2O[1][0];
        }
        else if( energies[i] > _stpowH2O[1][energies.size() - 1] ) {
            stpw_H2O_i = _stpowH2O[1][energies.size() - 1];
        }
        else {
            // Linear interpolation
            stpw_H2O_i = interp( energies[i] );
        }

        stopping_power[i] = stpw_H2O_i * density;
    }
pia.stammer's avatar
pia.stammer committed
204
    return stopping_power;
205
206
}

207
VectorVector Physics::GetTransportXS( Vector energies, Vector density ) {
208
    VectorVector transport_XS;
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
    double xsH2O_i;

    Spline interp;
    interp.set_points( _xsTransportH2O[0], _xsTransportH2O[1], false );    // false == linear interpolation

    for( unsigned i = 0; i < energies.size(); i++ ) {

        if( energies[i] < _xsTransportH2O[1][0] ) {
            xsH2O_i = _xsTransportH2O[1][0];
        }
        else if( energies[i] > _xsTransportH2O[1][energies.size() - 1] ) {
            xsH2O_i = _xsTransportH2O[1][energies.size() - 1];
        }
        else {
            // Linear interpolation
            xsH2O_i = interp( energies[i] );
        }

        transport_XS[i] = xsH2O_i * density;
    }
229
230
231
232
    return transport_XS;
}

Physics* Physics::Create() { return new Physics(); }
233

jannick.wolters's avatar
jannick.wolters committed
234
235
236
std::tuple<std::vector<VectorVector>, std::vector<VectorVector>> Physics::ReadENDL( std::string filename ) {
    std::vector<VectorVector> _headers;
    std::vector<VectorVector> _data;
steffen.schotthoefer's avatar
steffen.schotthoefer committed
237
238

    // open txt file
jannick.wolters's avatar
jannick.wolters committed
239
240
241
    std::string text_line, option_name;
    std::ifstream case_file;
    std::vector<std::string> option_value;
steffen.schotthoefer's avatar
steffen.schotthoefer committed
242
243
244

    /*--- Read the configuration file ---*/

jannick.wolters's avatar
jannick.wolters committed
245
    case_file.open( filename, std::ios::in );
steffen.schotthoefer's avatar
steffen.schotthoefer committed
246
247

    if( case_file.fail() ) {
248
        ErrorMessages::Error( "The ENDL file is missing!!", CURRENT_FUNCTION );
steffen.schotthoefer's avatar
steffen.schotthoefer committed
249
250
    }

jannick.wolters's avatar
jannick.wolters committed
251
    std::map<std::string, bool> included_options;
steffen.schotthoefer's avatar
steffen.schotthoefer committed
252
253
254
255

    /*--- Parse the physics file and save the values ---*/

    // list of entries (dataLine) of a datapackage, later becomes one entry in the list _data;
jannick.wolters's avatar
jannick.wolters committed
256
    std::list<Vector> dataPackList;
257
    VectorVector headerPack( 2 );
steffen.schotthoefer's avatar
steffen.schotthoefer committed
258
259
260
261
262
263
264
265

    // indices
    unsigned pack_idx        = 0;
    unsigned curr_HeaderLine = 0;

    while( getline( case_file, text_line ) ) {

        std::stringstream ss( text_line );    // give line to stringstream
jannick.wolters's avatar
jannick.wolters committed
266
        std::list<double> linelist;
steffen.schotthoefer's avatar
steffen.schotthoefer committed
267
268
269
270
271
        for( double double_in; ss >> double_in; ) {    // parse line
            linelist.push_back( double_in );
        }

        // Save content of line
272
        Vector dataLine( linelist.size() );
steffen.schotthoefer's avatar
steffen.schotthoefer committed
273
274
        bool header         = false;
        unsigned lineLenght = linelist.size();
steffen.schotthoefer's avatar
steffen.schotthoefer committed
275

steffen.schotthoefer's avatar
steffen.schotthoefer committed
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
        // check, if line is a dataline
        if( lineLenght == 2 || lineLenght == 3 ) {
            dataLine[0] = linelist.front();
            linelist.pop_front();
            dataLine[1] = linelist.front();
            linelist.pop_front();
        }
        if( lineLenght == 3 ) {
            dataLine[2] = linelist.front();
            linelist.pop_front();
        }
        // check, if line is  a header
        if( lineLenght == 8 || lineLenght == 9 ) {
            for( unsigned data_idx = 0; data_idx < lineLenght; data_idx++ ) {
                dataLine[data_idx] = linelist.front();
                linelist.pop_front();
                header = true;
            }
        }
        // check if line is an endline, then reset indices and save packs
        if( lineLenght == 1 ) {

            _headers.push_back( headerPack );

300
            VectorVector dataPack( dataPackList.size() );
steffen.schotthoefer's avatar
steffen.schotthoefer committed
301
302
303
304
305
306
307
308
309
            for( unsigned idx = 0; idx < dataPack.size(); idx++ ) {
                dataPack[idx] = dataPackList.front();
                dataPackList.pop_front();
            }
            _data.push_back( dataPack );

            pack_idx++;
            curr_HeaderLine = 0;
        }
steffen.schotthoefer's avatar
steffen.schotthoefer committed
310
311
312
313
314
315
316
317
318
319
        else {
            // save dataLine in vector;
            if( header ) {
                headerPack[curr_HeaderLine] = dataLine;
                curr_HeaderLine++;
            }
            else {
                dataPackList.push_back( dataLine );
            }
        }
steffen.schotthoefer's avatar
steffen.schotthoefer committed
320
    }
steffen.schotthoefer's avatar
steffen.schotthoefer committed
321
    case_file.close();
322
    return { _headers, _data };
steffen.schotthoefer's avatar
steffen.schotthoefer committed
323
}
324

325
326
VectorVector Physics::ReadStoppingPowers( std::string fileName ) {
    VectorVector stp_powers;
jannick.wolters's avatar
jannick.wolters committed
327
    std::string text_line;
328
329
330

    std::vector<double> e, stp_power;

jannick.wolters's avatar
jannick.wolters committed
331
332
    std::ifstream file_stream;
    file_stream.open( fileName, std::ios::in );
333
334
335

    if( file_stream.fail() ) {
        ErrorMessages::Error( "The Stopping Power file is missing!!", CURRENT_FUNCTION );
336
    }
337
338
339
340

    while( getline( file_stream, text_line ) ) {

        std::stringstream ss( text_line );    // give line to stringstream
jannick.wolters's avatar
jannick.wolters committed
341
        std::list<double> linelist;
342
343
344
        for( double double_in; ss >> double_in; ) {    // parse line
            linelist.push_back( double_in );
        }
345
346
347
348
349

        e.push_back( linelist.front() );    // append elements from line to column vectors
        linelist.pop_front();
        stp_power.push_back( linelist.front() );
        linelist.pop_front();
350
    }
pia.stammer's avatar
pia.stammer committed
351
    file_stream.close();
352

353
354
355
356
    stp_powers.push_back( Vector( e.size(), e.data() ) );
    stp_powers.push_back( Vector( stp_power.size(), stp_power.data() ) );

    return stp_powers;
steffen.schotthoefer's avatar
steffen.schotthoefer committed
357
}