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
00029 #ifndef G4NeutronHPFastLegendre_h
00030 #define G4NeutronHPFastLegendre_h 1
00031
00032 #include "globals.hh"
00033
00034 class G4NeutronHPFastLegendre
00035 {
00036 public:
00037
00038 G4NeutronHPFastLegendre()
00039 {
00040 value = new G4double * [31];
00041 value[0] = l0;
00042 value[1] = l1;
00043 value[2] = l2;
00044 value[3] = l3;
00045 value[4] = l4;
00046 value[5] = l5;
00047 value[6] = l6;
00048 value[7] = l7;
00049 value[8] = l8;
00050 value[9] = l9;
00051 value[10] = l10;
00052 value[11] = l11;
00053 value[12] = l12;
00054 value[13] = l13;
00055 value[14] = l14;
00056 value[15] = l15;
00057 value[16] = l16;
00058 value[17] = l17;
00059 value[18] = l18;
00060 value[19] = l19;
00061 value[20] = l20;
00062 value[21] = l21;
00063 value[22] = l22;
00064 value[23] = l23;
00065 value[24] = l24;
00066 value[25] = l25;
00067 value[26] = l26;
00068 value[27] = l27;
00069 value[28] = l28;
00070 value[29] = l29;
00071 value[30] = l30;
00072 integral = new G4double * [31];
00073 integral[0] = i0;
00074 integral[1] = i1;
00075 integral[2] = i2;
00076 integral[3] = i3;
00077 integral[4] = i4;
00078 integral[5] = i5;
00079 integral[6] = i6;
00080 integral[7] = i7;
00081 integral[8] = i8;
00082 integral[9] = i9;
00083 integral[10] = i10;
00084 integral[11] = i11;
00085 integral[12] = i12;
00086 integral[13] = i13;
00087 integral[14] = i14;
00088 integral[15] = i15;
00089 integral[16] = i16;
00090 integral[17] = i17;
00091 integral[18] = i18;
00092 integral[19] = i19;
00093 integral[20] = i20;
00094 integral[21] = i21;
00095 integral[22] = i22;
00096 integral[23] = i23;
00097 integral[24] = i24;
00098 integral[25] = i25;
00099 integral[26] = i26;
00100 integral[27] = i27;
00101 integral[28] = i28;
00102 integral[29] = i29;
00103 integral[30] = i30;
00104
00105 G4int i;
00106 for(i=0;i<31;i++) theNbin[i]=1+200*(i+1);
00107 }
00108
00109 ~G4NeutronHPFastLegendre()
00110 {
00111 delete [] value;
00112 delete [] integral;
00113 }
00114
00115 G4double Integrate(G4int l, G4double costh)
00116 {
00117 if(l>30) return regularIntegrate(l,costh);
00118 G4int bin = GetBin(l, costh);
00119 G4double y1, y2;
00120
00121 y1 = integral[l][bin];
00122 y2 = integral[l][bin+1];
00123
00124 return Interpolate(bin, l, y1, y2, costh);
00125 }
00126
00127 inline G4double Evaluate(G4int l, G4double costh)
00128 {
00129 if(l>30) return regularEvaluate(l,costh);
00130 G4double result;
00131 G4int bin = GetBin(l, costh);
00132 if(bin != theNbin[l]-1)
00133 {
00134 G4double y1, y2;
00135 y1 = value[l][bin];
00136 y2 = value[l][bin+1];
00137 result = Interpolate(bin, l, y1, y2, costh);
00138 }
00139 else
00140 {
00141 result = value[l][bin];
00142 }
00143 return result;
00144 }
00145
00146 private:
00147
00148 G4double regularEvaluate( int l , double x );
00149 G4double regularIntegrate( int l , double x );
00150
00151
00152 inline G4int GetBin(G4int l, G4double costh)
00153 {
00154 G4int bin=0;
00155 bin = G4int( (theNbin[l]-1)*(costh+1)/2. );
00156 if(bin == theNbin[l]-1) bin--;
00157 return bin;
00158 }
00159
00160 inline G4double Interpolate(G4int bin, G4int l, G4double y1, G4double y2, G4double x)
00161 {
00162 G4double slope = 0, off = 0, x2=0, x1mx2;
00163 G4int half = (theNbin[l]-1)/2;
00164
00165 x2 = (bin+1-half)/G4double(half);
00166 x1mx2 = 1./G4double( (theNbin[l]-1)/2 );
00167
00168 slope = (y2-y1)/x1mx2;
00169 off = y2-x2*slope;
00170 return x*slope+off;
00171 }
00172
00173 G4double ** value;
00174 G4double ** integral;
00175 G4int theNbin[31];
00176 static G4double l0[201];
00177 static G4double i0[201];
00178 static G4double l1[401];
00179 static G4double i1[401];
00180 static G4double l2[601];
00181 static G4double i2[601];
00182 static G4double l3[801];
00183 static G4double i3[801];
00184 static G4double l4[1001];
00185 static G4double i4[1001];
00186 static G4double l5[1201];
00187 static G4double i5[1201];
00188 static G4double l6[1401];
00189 static G4double i6[1401];
00190 static G4double l7[1601];
00191 static G4double i7[1601];
00192 static G4double l8[1801];
00193 static G4double i8[1801];
00194 static G4double l9[2001];
00195 static G4double i9[2001];
00196 static G4double l10[2201];
00197 static G4double i10[2201];
00198 static G4double l11[2401];
00199 static G4double i11[2401];
00200 static G4double l12[2601];
00201 static G4double i12[2601];
00202 static G4double l13[2801];
00203 static G4double i13[2801];
00204 static G4double l14[3001];
00205 static G4double i14[3001];
00206 static G4double l15[3201];
00207 static G4double i15[3201];
00208 static G4double l16[3401];
00209 static G4double i16[3401];
00210 static G4double l17[3601];
00211 static G4double i17[3601];
00212 static G4double l18[3801];
00213 static G4double i18[3801];
00214 static G4double l19[4001];
00215 static G4double i19[4001];
00216 static G4double l20[4201];
00217 static G4double i20[4201];
00218 static G4double l21[4401];
00219 static G4double i21[4401];
00220 static G4double l22[4601];
00221 static G4double i22[4601];
00222 static G4double l23[4801];
00223 static G4double i23[4801];
00224 static G4double l24[5001];
00225 static G4double i24[5001];
00226 static G4double l25[5201];
00227 static G4double i25[5201];
00228 static G4double l26[5401];
00229 static G4double i26[5401];
00230 static G4double l27[5601];
00231 static G4double i27[5601];
00232 static G4double l28[5801];
00233 static G4double i28[5801];
00234 static G4double l29[6001];
00235 static G4double i29[6001];
00236 static G4double l30[6201];
00237 static G4double i30[6201];
00238 };
00239 #endif