Geant4 Cross Reference |
1 // 2 // ******************************************************************** 3 // * License and Disclaimer * 4 // * * 5 // * The Geant4 software is copyright of the Copyright Holders of * 6 // * the Geant4 Collaboration. It is provided under the terms and * 7 // * conditions of the Geant4 Software License, included in the file * 8 // * LICENSE and available at http://cern.ch/geant4/license . These * 9 // * include a list of copyright holders. * 10 // * * 11 // * Neither the authors of this software system, nor their employing * 12 // * institutes,nor the agencies providing financial support for this * 13 // * work make any representation or warranty, express or implied, * 14 // * regarding this software system or assume any liability for its * 15 // * use. Please see the license in the file LICENSE and URL above * 16 // * for the full disclaimer and the limitation of liability. * 17 // * * 18 // * This code implementation is the result of the scientific and * 19 // * technical work of the GEANT4 collaboration. * 20 // * By using, copying, modifying or distributing the software (or * 21 // * any work based on the software) you agree to acknowledge its * 22 // * use in resulting scientific publications, and indicate your * 23 // * acceptance of all terms of the Geant4 Software license. * 24 // ******************************************************************** 25 // 26 27 #include "G4AdjointInterpolator.hh" 28 29 G4ThreadLocal G4AdjointInterpolator* G4AdjointInterpolator::fInstance = nullptr; 30 31 /////////////////////////////////////////////////////// 32 G4AdjointInterpolator* G4AdjointInterpolator::GetAdjointInterpolator() 33 { 34 return GetInstance(); 35 } 36 37 /////////////////////////////////////////////////////// 38 G4AdjointInterpolator* G4AdjointInterpolator::GetInstance() 39 { 40 if(!fInstance) 41 { 42 fInstance = new G4AdjointInterpolator; 43 } 44 return fInstance; 45 } 46 47 /////////////////////////////////////////////////////// 48 G4AdjointInterpolator::G4AdjointInterpolator() {} 49 50 /////////////////////////////////////////////////////// 51 G4AdjointInterpolator::~G4AdjointInterpolator() {} 52 53 /////////////////////////////////////////////////////// 54 G4double G4AdjointInterpolator::LinearInterpolation(G4double& x, G4double& x1, 55 G4double& x2, G4double& y1, 56 G4double& y2) 57 { 58 G4double res = y1 + (x - x1) * (y2 - y1) / (x2 - x1); 59 return res; 60 } 61 62 /////////////////////////////////////////////////////// 63 G4double G4AdjointInterpolator::LogarithmicInterpolation( 64 G4double& x, G4double& x1, G4double& x2, G4double& y1, G4double& y2) 65 { 66 if(y1 <= 0. || y2 <= 0. || x1 <= 0.) 67 return LinearInterpolation(x, x1, x2, y1, y2); 68 G4double B = std::log(y2 / y1) / std::log(x2 / x1); 69 G4double A = y1 / std::pow(x1, B); 70 G4double res = A * std::pow(x, B); 71 return res; 72 } 73 74 /////////////////////////////////////////////////////// 75 G4double G4AdjointInterpolator::ExponentialInterpolation( 76 G4double& x, G4double& x1, G4double& x2, G4double& y1, G4double& y2) 77 { 78 G4double B = (std::log(y2) - std::log(y1)) / (x2 - x1); 79 G4double A = y1 * std::exp(-B * x1); 80 G4double res = A * std::exp(B * x); 81 return res; 82 } 83 84 /////////////////////////////////////////////////////// 85 G4double G4AdjointInterpolator::Interpolation(G4double& x, G4double& x1, 86 G4double& x2, G4double& y1, 87 G4double& y2, 88 const G4String& InterPolMethod) 89 { 90 if(InterPolMethod == "Log") 91 { 92 return LogarithmicInterpolation(x, x1, x2, y1, y2); 93 } 94 else if(InterPolMethod == "Lin") 95 { 96 return LinearInterpolation(x, x1, x2, y1, y2); 97 } 98 else if(InterPolMethod == "Exp") 99 { 100 return ExponentialInterpolation(x, x1, x2, y1, y2); 101 } 102 else 103 { 104 G4ExceptionDescription ed; 105 ed << "The interpolation method that you invoked does not exist!\n"; 106 G4Exception("G4AdjointInterpolator::Interpolation", "adoint001", 107 FatalException, ed); 108 return 0.; 109 } 110 } 111 112 /////////////////////////////////////////////////////// 113 // only valid if x_vec is monotically increasing 114 std::size_t G4AdjointInterpolator::FindPosition(G4double& x, 115 std::vector<G4double>& x_vec, std::size_t, 116 std::size_t) 117 { 118 // most rapid method could be used probably 119 120 std::size_t ndim = x_vec.size(); 121 std::size_t ind1 = 0; 122 std::size_t ind2 = ndim - 1; 123 124 if(ndim > 1) 125 { 126 if(x_vec[0] < x_vec[1]) 127 { // increasing 128 do 129 { 130 std::size_t midBin = (ind1 + ind2) / 2; 131 if(x < x_vec[midBin]) 132 ind2 = midBin; 133 else 134 ind1 = midBin; 135 } while(ind2 - ind1 > 1); 136 } 137 else 138 { 139 do 140 { 141 std::size_t midBin = (ind1 + ind2) / 2; 142 if(x < x_vec[midBin]) 143 ind1 = midBin; 144 else 145 ind2 = midBin; 146 } while(ind2 - ind1 > 1); 147 } 148 } 149 150 return ind1; 151 } 152 153 /////////////////////////////////////////////////////// 154 // only valid if x_vec is monotically increasing 155 std::size_t G4AdjointInterpolator::FindPositionForLogVector( 156 G4double& log_x, std::vector<G4double>& log_x_vec) 157 { 158 // most rapid method could be used probably 159 return FindPosition(log_x, log_x_vec); 160 } 161 162 /////////////////////////////////////////////////////// 163 G4double G4AdjointInterpolator::Interpolate(G4double& x, 164 std::vector<G4double>& x_vec, 165 std::vector<G4double>& y_vec, 166 const G4String& InterPolMethod) 167 { 168 std::size_t i = FindPosition(x, x_vec); 169 return Interpolation(x, x_vec[i], x_vec[i + 1], y_vec[i], y_vec[i + 1], 170 InterPolMethod); 171 } 172 173 /////////////////////////////////////////////////////// 174 G4double G4AdjointInterpolator::InterpolateWithIndexVector( 175 G4double& x, std::vector<G4double>& x_vec, std::vector<G4double>& y_vec, 176 std::vector<std::size_t>& index_vec, G4double x0, 177 G4double dx) // only linear interpolation possible 178 { 179 std::size_t ind = 0; 180 if(x > x0) 181 ind = G4int((x - x0) / dx); 182 if(ind >= index_vec.size() - 1) 183 ind = index_vec.size() - 2; 184 std::size_t ind1 = index_vec[ind]; 185 std::size_t ind2 = index_vec[ind + 1]; 186 if(ind1 > ind2) 187 { 188 std::size_t ind11 = ind1; 189 ind1 = ind2; 190 ind2 = ind11; 191 } 192 ind = FindPosition(x, x_vec, ind1, ind2); 193 return Interpolation(x, x_vec[ind], x_vec[ind + 1], y_vec[ind], 194 y_vec[ind + 1], "Lin"); 195 } 196 197 /////////////////////////////////////////////////////// 198 G4double G4AdjointInterpolator::InterpolateForLogVector( 199 G4double& log_x, std::vector<G4double>& log_x_vec, 200 std::vector<G4double>& log_y_vec) 201 { 202 std::size_t i = FindPositionForLogVector(log_x, log_x_vec); 203 204 G4double log_y = LinearInterpolation(log_x, log_x_vec[i], log_x_vec[i + 1], 205 log_y_vec[i], log_y_vec[i + 1]); 206 return log_y; 207 } 208