rtang/ADPlusPlus
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity

simplify Qk.h

Browse Source
This commit is contained in:
Ryan@iMac 2022-10-12 15:26:17 -04:00
parent 78576a70d2
commit ddfbf67404
1 changed files with 22 additions and 40 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

62
Qk.h
View File

@ -40,61 +40,43 @@ double * QK(double Energy_keV, double radius_cm, double distance_cm, double thic
double delx2 = (gamma-alpha)/div;
for(int i = 0; i<=div; i++){
int J = i % 2;
if( J == 0 ){
A = 2.0;
}else{
A = 4.0;
beta1 = i * delx1;
beta2 = alpha + i * delx2;
}
if( i == 0 || i == div ) {
A = 1.0;
beta1 = i * delx1;
beta2 = alpha + i * delx2;
}else{
A = ( i % 2 == 0 ? 2.0 : 4.0 );
if( i % 2 == 1 ){
beta1 = i * delx1;
beta2 = alpha + i * delx2;
}
}
double cosb = cos(beta1);
double sinb = sin(beta1);
double ex1 = exp( - 1.0 * Tau * thickness_cm / cosb );
double term3 = (1-ex1) * sinb * A * delx1;
double term1 = LegendreP(2, beta1) * term3;
double term2 = LegendreP(4, beta1) * term3;
double ex1 = exp( - 1.0 * Tau * thickness_cm / cos(beta1) );
double term1 = A * (1-ex1) * sin(beta1) * delx1;
sum1 += term1;
sum2 += term2;
sum3 += term3;
sum1 += LegendreP(2, beta1) * term1;
sum2 += LegendreP(4, beta1) * term1;
sum3 += term1;
cosb = cos(beta2);
sinb = sin(beta2);
double ex2 = exp( -1 * Tau * (radius_cm / sinb - distance_cm /cosb) );
double term6 = A * (1-ex2) * sinb * delx2;
double term4 = LegendreP(2, beta2) * term6;
double term5 = LegendreP(4, beta2) * term6;
double ex2 = exp( -1 * Tau * (radius_cm / sin(beta2) - distance_cm /cos(beta2)) );
double term2 = A * (1-ex2) * sin(beta2) * delx2;
sum4 = sum4 +term4;
sum5 = sum5 +term5;
sum6 = sum6 +term6;
sum4 += LegendreP(2, beta2) * term2;
sum5 += LegendreP(4, beta2) * term2;
sum6 += term2;
//if( i % 75 == 0) printf("%4d | %10.6f, %10.6f %d | %10.6f, %10.6f %10.6f \n", i, beta, A, J, sum1, sum2, sum3);
//if( i % 75 == 0) printf("%4d | %10.6f, %10.6f, %10.6f %d | %10.6f, %10.6f %10.6f %10.6f, %10.6f %10.6f \n", i, beta1, beta2, A, J, sum1, sum2, sum3, sum4, sum5, sum6);
}
double ans1 = sum1/3;
double ans2 = sum2/3;
double ans3 = sum3/3;
double ans4 = sum4/3;
double ans5 = sum5/3;
double ans6 = sum6/3;
double * Qk = new double[2];
Qk[0] = (ans1+ans4)/(ans3+ans6);
Qk[1] = (ans2+ans5)/(ans3+ans6);
Qk[0] = (sum1+sum4)/(sum3+sum6);
Qk[1] = (sum2+sum5)/(sum3+sum6);
printf("--------------\n");
printf(" QD2 = %lf\n",Qk[0]);
printf(" QD4 = %lf\n",Qk[1]);
printf(" QD2 = %.20f\n",Qk[0]);
printf(" QD4 = %.20f\n",Qk[1]);
printf("--------------\n");
return Qk;