00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028 #include "G4AdjointCSMatrix.hh"
00029 #include "G4AdjointInterpolator.hh"
00030
00031 G4AdjointInterpolator* G4AdjointInterpolator::theInstance = 0;
00033
00034 G4AdjointInterpolator* G4AdjointInterpolator::GetAdjointInterpolator()
00035 { if(theInstance == 0) {
00036 static G4AdjointInterpolator interpolator;
00037 theInstance = &interpolator;
00038 }
00039 return theInstance;
00040 }
00042
00043 G4AdjointInterpolator* G4AdjointInterpolator::GetInstance()
00044 { if(theInstance == 0) {
00045 static G4AdjointInterpolator interpolator;
00046 theInstance = &interpolator;
00047 }
00048 return theInstance;
00049 }
00050
00052
00053 G4AdjointInterpolator::G4AdjointInterpolator()
00054 {;
00055 }
00057
00058 G4AdjointInterpolator::~G4AdjointInterpolator()
00059 {;
00060 }
00062
00063 G4double G4AdjointInterpolator::LinearInterpolation(G4double& x,G4double& x1,G4double& x2,G4double& y1,G4double& y2)
00064 { G4double res = y1+ (x-x1)*(y2-y1)/(x2-x1);
00065
00066 return res;
00067 }
00069
00070 G4double G4AdjointInterpolator::LogarithmicInterpolation(G4double& x,G4double& x1,G4double& x2,G4double& y1,G4double& y2)
00071 { if (y1<=0 || y2<=0 || x1<=0) return LinearInterpolation(x,x1,x2,y1,y2);
00072 G4double B=std::log(y2/y1)/std::log(x2/x1);
00073
00074 G4double A=y1/std::pow(x1,B);
00075 G4double res=A*std::pow(x,B);
00076
00077 return res;
00078 }
00080
00081 G4double G4AdjointInterpolator::ExponentialInterpolation(G4double& x,G4double& x1,G4double& x2,G4double& y1,G4double& y2)
00082 { G4double B=(std::log(y2)-std::log(y1));
00083 B=B/(x2-x1);
00084 G4double A=y1*std::exp(-B*x1);
00085 G4double res=A*std::exp(B*x);
00086 return res;
00087 }
00089
00090 G4double G4AdjointInterpolator::Interpolation(G4double& x,G4double& x1,G4double& x2,G4double& y1,G4double& y2,G4String InterPolMethod)
00091 {
00092 if (InterPolMethod == "Log" ){
00093 return LogarithmicInterpolation(x,x1,x2,y1,y2);
00094 }
00095 else if (InterPolMethod == "Lin" ){
00096 return LinearInterpolation(x,x1,x2,y1,y2);
00097 }
00098 else if (InterPolMethod == "Exp" ){
00099 return ExponentialInterpolation(x,x1,x2,y1,y2);
00100 }
00101 else {
00102
00103 return -1111111111.;
00104 }
00105 }
00107
00108 size_t G4AdjointInterpolator::FindPosition(G4double& x,std::vector<G4double>& x_vec,size_t , size_t )
00109 {
00110
00111
00112
00113 size_t ndim = x_vec.size();
00114 size_t ind1 = 0;
00115 size_t ind2 = ndim - 1;
00116
00117
00118
00119
00120
00121
00122
00123
00124
00125 if (ndim >1) {
00126
00127 if (x_vec[0] < x_vec[1] ) {
00128 do {
00129 size_t midBin = (ind1 + ind2)/2;
00130 if (x < x_vec[midBin])
00131 ind2 = midBin;
00132 else
00133 ind1 = midBin;
00134 } while (ind2 - ind1 > 1);
00135 }
00136 else {
00137 do {
00138 size_t midBin = (ind1 + ind2)/2;
00139 if (x < x_vec[midBin])
00140 ind1 = midBin;
00141 else
00142 ind2 = midBin;
00143 } while (ind2 - ind1 > 1);
00144 }
00145
00146 }
00147
00148 return ind1;
00149 }
00150
00152
00153 size_t G4AdjointInterpolator::FindPositionForLogVector(G4double& log_x,std::vector<G4double>& log_x_vec)
00154 {
00155
00156 return FindPosition(log_x, log_x_vec);
00157
00158
00159
00160
00161
00162
00163
00164
00165
00166
00167
00168
00169
00170
00171
00172
00173 }
00175
00176 G4double G4AdjointInterpolator::Interpolate(G4double& x,std::vector<G4double>& x_vec,std::vector<G4double>& y_vec,G4String InterPolMethod)
00177 { size_t i=FindPosition(x,x_vec);
00178
00179
00180
00181 return Interpolation( x,x_vec[i],x_vec[i+1],y_vec[i],y_vec[i+1],InterPolMethod);
00182 }
00183
00185
00186 G4double G4AdjointInterpolator::InterpolateWithIndexVector(G4double& x,std::vector<G4double>& x_vec,std::vector<G4double>& y_vec,
00187 std::vector<size_t>& index_vec,G4double x0, G4double dx)
00188 { size_t ind=0;
00189 if (x>x0) ind=int((x-x0)/dx);
00190 if (ind >= index_vec.size()-1) ind= index_vec.size()-2;
00191 size_t ind1 = index_vec[ind];
00192 size_t ind2 = index_vec[ind+1];
00193 if (ind1 >ind2) {
00194 size_t ind11=ind1;
00195 ind1=ind2;
00196 ind2=ind11;
00197
00198 }
00199
00200 ind=FindPosition(x,x_vec,ind1,ind2);
00201 return Interpolation( x,x_vec[ind],x_vec[ind+1],y_vec[ind],y_vec[ind+1],"Lin");
00202
00203 }
00204
00205
00207
00208 G4double G4AdjointInterpolator::InterpolateForLogVector(G4double& log_x,std::vector<G4double>& log_x_vec,std::vector<G4double>& log_y_vec)
00209 {
00210 size_t i=FindPositionForLogVector(log_x,log_x_vec);
00211
00212
00213
00214
00215
00216
00217
00218
00219 G4double log_y=LinearInterpolation(log_x,log_x_vec[i],log_x_vec[i+1],log_y_vec[i],log_y_vec[i+1]);
00220 return log_y;
00221
00222 }