#include math.h
#define PI 3.14159265358979323846
#define abs(x) ((x)>0.?(x):-(x))
#define sgn(x) ((x)>0.?1.:-1.)
#define NSUB 10
#define twoPI (3.14159265358979323846*2.0)
double zf,r,ya,yb,za,zb;
double bigN;
double ra,rb;
main(argc,argv)
int argc;
char *argv[];
{
double Ir(),IrM();
double z,Rab,e1s,e2s,e1o,e2o,Ha,Hb,ha,hb,theta,dy,ES,ESM;
double rA, rB, ES0;
int i,j,loc,nloc,NZ,Ntheta;
char buf[80];
FILE *fp;
if (argc != 10){
fprintf(stderr,"Usage: %s ra rb Rab Ntheta e1s e2s e1o e2o NZ\n", argv[0]);
exit(1);
}
rA = atof(argv[1]);
rB = atof(argv[2]);
Rab = atof(argv[3]);
Ntheta = atoi(argv[4]);
e1s = atof(argv[5]);
e2s = atof(argv[6]);
e1o = atof(argv[7]);
e2o = atof(argv[8]);
NZ = atoi(argv[9]);
ES0 = 0.5*(1.0/rA+1.0/rB)*(1.0/e1o-1.0/e1s);
if (Rab < rA+rB){
fprintf(stderr,"Rab = %f must be greater than rA+rB = %f\n",Rab,rA+rB);
exit(1);
}
for (j=0;j<=Ntheta;j++){
printf("%-7.2f",j*90/(Ntheta-1.0));
fprintf(stderr,"%-7.2f",j*90/(Ntheta*1.0));
}
printf("\n");
fprintf(stderr,"\n");
/*Main loop over the location of the center point of the ion-pair*/ for (i=0;i<=NZ;i++){
z = -8+16.0*(i+0.333)/(1.0*NZ);
printf("%-7.2f",z);
fprintf(stderr,"%-7.2f",z);
for (j=0;j<=Ntheta;j++){
theta = j*PI/(2*(Ntheta-1.0));
dy = Rab*sin(theta);
Ha = z-0.5*Rab*cos(theta);
Hb = z+0.5*Rab*cos(theta);
nloc = 0;
if (Ha <= -rA && Hb <= -rB) {
loc = 1;
nloc++;}
if (Ha <= -rA && Hb <= 0 && Hb > -rB) {
loc = 2;
nloc++;}
if (Ha <= -rA && Hb <= rB && Hb > 0) {
loc = 3;
nloc++;}
if (Ha <= -rA && Hb > rB) {
loc = 4;
nloc++;}
if (Ha <= 0 && Ha > -rA && Hb > rB) {
loc = 5;
nloc++;}
if (Ha <= rA && Ha > 0 && Hb > rB) {
loc = 6;
nloc++;}
if (Ha > rA && Hb > rB) {
loc = 7;
nloc++;}
if (Ha <= 0 && Ha > -rA && Hb <= 0 && Hb > -rB) {
loc = 8;
nloc++;}
if (Ha <= 0 && Ha > -rA && Hb <= rB && Hb > 0) {
loc = 9;
nloc++;}
if (Ha <= rA && Ha > 0 && Hb <= rB && Hb >0 ) {
loc = 10;
nloc++;}
if (nloc == 0){
printf("location undefined Ha = %f Hb = %f\n",Ha,Hb);
}
if (nloc > 1){
printf("location defined %d times. Ha = %f, Hb = %f\n",nloc,Ha,Hb);
}
if (loc < 5){
ra = rA;
rb = rB;
ha = abs(Ha);
hb = abs(Hb);
ES = Ir(e1o,e2o,ha,hb,dy,loc)-Ir(e1s,e2s,ha,hb,dy,loc);
if (loc == 4){
/*Marcus*/ ESM = IrM(e1o,e2o,ha,hb,Rab)-IrM(e1s,e2s,ha,hb,Rab);
fprintf(stderr,"%-10.6f",ESM/ES0);
}
}
if (loc > 4 && loc < 8){
ra = rB;
rb = rA;
ha = abs(Hb);
hb = abs(Ha);
ES = Ir(e2o,e1o,ha,hb,dy,8-loc)-Ir(e2s,e1s,ha,hb,dy,8-loc);
}
if (loc >7 && loc <10){
ra = rA;
rb = rB;
ha = abs(Ha);
hb = abs(Hb);
ES = Ir(e1o,e2o,ha,hb,dy,loc-3)-Ir(e1s,e2s,ha,hb,dy,loc-3);
}
if (loc == 10){
ra = rB;
rb = rA;
ha = abs(Hb);
hb = abs(Ha);
ES = Ir(e2o,e1o,ha,hb,dy,15-loc)-Ir(e2s,e1s,ha,hb,dy,15-loc);
}
printf("%-10.6f",ES/ES0);
}
printf("\n");
fprintf(stderr,"\n");
}
}
double IrM(e1,e2,ha,hb,R) double e1,e2,ha,hb,R;
{
double eta,val;
eta = (e1-e2)/(e1+e2);
val = 0.5/(ra*e1)+0.5/(rb*e2)+0.25*eta/(e1*ha)-0.25*eta/(e2*hb)-2.0/(R*(e1+e2));
return (val);
}
double Ir(e1,e2,ha,hb,dy,n) double e1,e2,ha,hb,dy;
int n;
{
double eta,val,valA,valB, valA1, valA2, valB1,valB2,dsint();
eta = (e1-e2)/(e1+e2);
switch (n){
case 1: /*both ions are totaly in the e1 region (z<0)*/
valA = dsint(-ra,ra,0.,0.,0.,0.,ra) + eta*dsint(-ra,ra,2*ha,0.,0.,0.,ra) - dsint(-ra,ra,ha-hb,0.,-dy,0.,ra) - eta*dsint(-ra,ra,ha+hb,0.,-dy,0.,ra) + eta*dsint(-ra,ra,0.,2*ha,0.,0.,ra) + eta*eta*dsint(-ra,ra,2*ha,2*ha,0.,0.,ra) - eta*dsint(-ra,ra,ha-hb,2*ha,-dy,0.,ra) - eta*eta*dsint(-ra,ra,ha+hb,2*ha,-dy,0.,ra) - dsint(-ra,ra,0.,ha-hb,0.,-dy,ra) - eta*dsint(-ra,ra,2*ha,ha-hb,0.,-dy,ra) + dsint(-ra,ra,ha-hb,ha-hb,-dy,-dy,ra) + eta*dsint(-ra,ra,ha+hb,ha-hb,-dy,-dy,ra) - eta*dsint(-ra,ra,0.,ha+hb,0.,-dy,ra) - eta*eta*dsint(-ra,ra,2*ha,ha+hb,0.,-dy,ra) + eta*dsint(-ra,ra,ha-hb,ha+hb,-dy,-dy,ra) + eta*eta*dsint(-ra,ra,ha+hb,ha+hb,-dy,-dy,ra);
valA /= 8*PI*e1;
valB = dsint(-rb,rb,hb-ha,hb-ha,dy,dy,rb) + eta*dsint(-rb,rb,ha+hb,hb-ha,dy,dy,rb) - dsint(-rb,rb,0.,hb-ha,0.,dy,rb) - eta*dsint(-rb,rb,2*hb,hb-ha,0.,dy,rb) + eta*dsint(-rb,rb,hb-ha,ha+hb,dy,dy,rb) + eta*eta*dsint(-rb,rb,ha+hb,ha+hb,dy,dy,rb) - eta*dsint(-rb,rb,0.,ha+hb,0.,dy,rb) - eta*eta*dsint(-rb,rb,2*hb,ha+hb,0.,dy,rb) - dsint(-rb,rb,hb-ha,0.,dy,0.,rb) - eta*dsint(-rb,rb,ha+hb,0.,dy,0.,rb) + dsint(-rb,rb,0.,0.,0.,0.,rb) + eta*dsint(-rb,rb,2*hb,0.,0.,0.,rb) - eta*dsint(-rb,rb,hb-ha,2*hb,dy,0.,rb) - eta*eta*dsint(-rb,rb,ha+hb,2*hb,dy,0.,rb) + eta*dsint(-rb,rb,0.,2*hb,0.,0.,rb) + eta*eta*dsint(-rb,rb,2*hb,2*hb,0.,0.,rb);
valB /= 8*PI*e1;
val = valA + valB;
break;
case 2: /*both centers are in e1, one ion is in e1 2nd ion is partially in e1*/
valA = dsint(-ra,ra,0.,0.,0.,0.,ra) + eta*dsint(-ra,ra,2*ha,0.,0.,0.,ra) - dsint(-ra,ra,ha-hb,0.,-dy,0.,ra) - eta*dsint(-ra,ra,ha+hb,0.,-dy,0.,ra) + eta*dsint(-ra,ra,0.,2*ha,0.,0.,ra) + eta*eta*dsint(-ra,ra,2*ha,2*ha,0.,0.,ra) - eta*dsint(-ra,ra,ha-hb,2*ha,-dy,0.,ra) - eta*eta*dsint(-ra,ra,ha+hb,2*ha,-dy,0.,ra) - dsint(-ra,ra,0.,ha-hb,0.,-dy,ra) - eta*dsint(-ra,ra,2*ha,ha-hb,0.,-dy,ra) + dsint(-ra,ra,ha-hb,ha-hb,-dy,-dy,ra) + eta*dsint(-ra,ra,ha+hb,ha-hb,-dy,-dy,ra) - eta*dsint(-ra,ra,0.,ha+hb,0.,-dy,ra) - eta*eta*dsint(-ra,ra,2*ha,ha+hb,0.,-dy,ra) + eta*dsint(-ra,ra,ha-hb,ha+hb,-dy,-dy,ra) + eta*eta*dsint(-ra,ra,ha+hb,ha+hb,-dy,-dy,ra);
valA /= 8*PI*e1;
valB1 = dsint(-hb,rb,hb-ha,hb-ha,dy,dy,rb) + eta*dsint(-hb,rb,ha+hb,hb-ha,dy,dy,rb) - dsint(-hb,rb,0.,hb-ha,0.,dy,rb) - eta*dsint(-hb,rb,2*hb,hb-ha,0.,dy,rb) + eta*dsint(-hb,rb,hb-ha,ha+hb,dy,dy,rb) + eta*eta*dsint(-hb,rb,ha+hb,ha+hb,dy,dy,rb) - eta*dsint(-hb,rb,0.,ha+hb,0.,dy,rb) - eta*eta*dsint(-hb,rb,2*hb,ha+hb,0.,dy,rb) - dsint(-hb,rb,hb-ha,0.,dy,0.,rb) - eta*dsint(-hb,rb,ha+hb,0.,dy,0.,rb) + dsint(-hb,rb,0.,0.,0.,0.,rb) + eta*dsint(-hb,rb,2*hb,0.,0.,0.,rb) - eta*dsint(-hb,rb,hb-ha,2*hb,dy,0.,rb) - eta*eta*dsint(-hb,rb,ha+hb,2*hb,dy,0.,rb) + eta*dsint(-hb,rb,0.,2*hb,0.,0.,rb) + eta*eta*dsint(-hb,rb,2*hb,2*hb,0.,0.,rb);
valB1 /= 8*PI*e1;
valB2 = dsint(-rb,-hb,hb-ha,hb-ha,dy,dy,rb) - dsint(-rb,-hb,hb-ha,0.,dy,0.,rb) - dsint(-rb,-hb,0.,hb-ha,0.,dy,rb) + dsint(-rb,-hb,0.,0.,0.,0.,rb);
valB2 *= (e2/(2*PI*(e1+e2)*(e1+e2)));
val = valA + valB1 + valB2;
break;
case 3: /*one ion is in e1 2nd ion is partially in e1, but its center is in e2*/
valA = dsint(-ra,ra,0.,0.,0.,0.,ra) + eta*dsint(-ra,ra,2*ha,0.,0.,0.,ra) - (1+eta)*dsint(-ra,ra,ha+hb,0.,-dy,0.,ra) + eta*dsint(-ra,ra,0.,2*ha,0.,0.,ra) + eta*eta*dsint(-ra,ra,2*ha,2*ha,0.,0.,ra) - eta*(1+eta)*dsint(-ra,ra,ha+hb,2*ha,-dy,0.,ra) - (1+eta)*dsint(-ra,ra,0.,ha+hb,0.,-dy,ra) - eta*(1+eta)*dsint(-ra,ra,2*ha,ha+hb,0.,-dy,ra) + (1+eta)*(1+eta)*dsint(-ra,ra,ha+hb,ha+hb,-dy,-dy,ra);
valA /= 8*PI*e1;
valB1 = dsint(hb,rb,-ha-hb,-ha-hb,dy,dy,rb) + eta*dsint(hb,rb,ha-hb,-ha-hb,dy,dy,rb) - (1+eta)*dsint(hb,rb,0.,-ha-hb,0.,dy,rb) + eta*dsint(hb,rb,-ha-hb,ha-hb,dy,dy,rb) + eta*eta*dsint(hb,rb,ha-hb,ha-hb,dy,dy,rb) - eta*(1+eta)*dsint(hb,rb,0.,ha-hb,0.,dy,rb) - (1+eta)*dsint(hb,rb,-ha-hb,0.,dy,0.,rb) - eta*(1+eta)*dsint(hb,rb,ha-hb,0.,dy,0.,rb) + (1+eta)*(1+eta)*dsint(hb,rb,0.,0.,0.,0.,rb);
valB1 /= 8*PI*e1;
valB2 = dsint(-rb,hb,0.,0.,0.,0.,rb) - eta*dsint(-rb,hb,-2*hb,0.,0.,0.,rb) - (1-eta)*dsint(-rb,hb,-ha-hb,0.,dy,0.,rb) - eta*dsint(-rb,hb,0.,-2*hb,0.,0.,rb) + eta*eta*dsint(-rb,hb,-2*hb,-2*hb,0.,0.,rb) + eta*(1-eta)*dsint(-rb,hb,-ha-hb,-2*hb,dy,0.,rb) - (1-eta)*dsint(-rb,hb,0.,-ha-hb,0.,dy,rb) + eta*(1-eta)*dsint(-rb,hb,-2*hb,-ha-hb,0.,dy,rb) + (1-eta)*(1-eta)*dsint(-rb,hb,-ha-hb,-ha-hb,dy,dy,rb);
valB2 /= 8*PI*e2;
val = valA + valB1 + valB2;
break;
case 4: /*one ion is totaly in e1 2nd ion is totaly in e2*/
valA = dsint(-ra,ra,0.,0.,0.,0.,ra) + eta*dsint(-ra,ra,2*ha,0.,0.,0.,ra) - (1+eta)*dsint(-ra,ra,ha+hb,0.,-dy,0.,ra) + eta*dsint(-ra,ra,0.,2*ha,0.,0.,ra) + eta*eta*dsint(-ra,ra,2*ha,2*ha,0.,0.,ra) - eta*(1+eta)*dsint(-ra,ra,ha+hb,2*ha,-dy,0.,ra) - (1+eta)*dsint(-ra,ra,0.,ha+hb,0.,-dy,ra) - eta*(1+eta)*dsint(-ra,ra,2*ha,ha+hb,0.,-dy,ra) + (1+eta)*(1+eta)*dsint(-ra,ra,ha+hb,ha+hb,-dy,-dy,ra);
valA /= 8*PI*e1;
valB = dsint(-rb,rb,0.,0.,0.,0.,rb) - eta*dsint(-rb,rb,-2*hb,0.,0.,0.,rb) - (1-eta)*dsint(-rb,rb,-ha-hb,0.,dy,0.,rb) - eta*dsint(-rb,rb,0.,-2*hb,0.,0.,rb) + eta*eta*dsint(-rb,rb,-2*hb,-2*hb,0.,0.,rb) + eta*(1-eta)*dsint(-rb,rb,-ha-hb,-2*hb,dy,0.,rb) - (1-eta)*dsint(-rb,rb,0.,-ha-hb,0.,dy,rb) + eta*(1-eta)*dsint(-rb,rb,-2*hb,-ha-hb,0.,dy,rb) + (1-eta)*(1-eta)*dsint(-rb,rb,-ha-hb,-ha-hb,dy,dy,rb);
valB /= 8*PI*e2;
val = valA + valB;
break;
case 5: /*Both ions are partially in e1 with both centers in e1*/
valA1 = dsint(-ha,ra,0.,0.,0.,0.,ra) + eta*dsint(-ha,ra,2*ha,0.,0.,0.,ra) - dsint(-ha,ra,ha-hb,0.,-dy,0.,ra) - eta*dsint(-ha,ra,ha+hb,0.,-dy,0.,ra) + eta*dsint(-ha,ra,0.,2*ha,0.,0.,ra) + eta*eta*dsint(-ha,ra,2*ha,2*ha,0.,0.,ra) - eta*dsint(-ha,ra,ha-hb,2*ha,-dy,0.,ra) - eta*eta*dsint(-ha,ra,ha+hb,2*ha,-dy,0.,ra) - dsint(-ha,ra,0.,ha-hb,0.,-dy,ra) - eta*dsint(-ha,ra,2*ha,ha-hb,0.,-dy,ra) + dsint(-ha,ra,ha-hb,ha-hb,-dy,-dy,ra) + eta*dsint(-ha,ra,ha+hb,ha-hb,-dy,-dy,ra) - eta*dsint(-ha,ra,0.,ha+hb,0.,-dy,ra) - eta*eta*dsint(-ha,ra,2*ha,ha+hb,0.,-dy,ra) + eta*dsint(-ha,ra,ha-hb,ha+hb,-dy,-dy,ra) + eta*eta*dsint(-ha,ra,ha+hb,ha+hb,-dy,-dy,ra);
valA1 /= 8*PI*e1;
valA2 = dsint(-ra,-ha,0.,0.,0.,0.,ra) - dsint(-ra,-ha,ha-hb,0.,-dy,0.,ra) - dsint(-ra,-ha,0.,ha-hb,0.,-dy,ra) + dsint(-ra,-ha,ha-hb,ha-hb,-dy,-dy,ra);
valA2 *= (e2/(2*PI*(e1+e2)*(e1+e2)));
valB1 = dsint(-hb,rb,hb-ha,hb-ha,dy,dy,rb) + eta*dsint(-hb,rb,ha+hb,hb-ha,dy,dy,rb) - dsint(-hb,rb,0.,hb-ha,0.,dy,rb) - eta*dsint(-hb,rb,2*hb,hb-ha,0.,dy,rb) + eta*dsint(-hb,rb,hb-ha,ha+hb,dy,dy,rb) + eta*eta*dsint(-hb,rb,ha+hb,ha+hb,dy,dy,rb) - eta*dsint(-hb,rb,0.,ha+hb,0.,dy,rb) - eta*eta*dsint(-hb,rb,2*hb,ha+hb,0.,dy,rb) - dsint(-hb,rb,hb-ha,0.,dy,0.,rb) - eta*dsint(-hb,rb,ha+hb,0.,dy,0.,rb) + dsint(-hb,rb,0.,0.,0.,0.,rb) + eta*dsint(-hb,rb,2*hb,0.,0.,0.,rb) - eta*dsint(-hb,rb,hb-ha,2*hb,dy,0.,rb) - eta*eta*dsint(-hb,rb,ha+hb,2*hb,dy,0.,rb) + eta*dsint(-hb,rb,0.,2*hb,0.,0.,rb) + eta*eta*dsint(-hb,rb,2*hb,2*hb,0.,0.,rb);
valB1 /= 8*PI*e1;
valB2 = dsint(-rb,-hb,hb-ha,hb-ha,dy,dy,rb) - dsint(-rb,-hb,hb-ha,0.,dy,0.,rb) - dsint(-rb,-hb,0.,hb-ha,0.,dy,rb) + dsint(-rb,-hb,0.,0.,0.,0.,rb);
valB2 *= (e2/(2*PI*(e1+e2)*(e1+e2)));
val = valA1 + valA2 + valB1 + valB2;
break;
case 6: /*Both ions are partially in e1 one center in e1, the second in e2*/
valA1 = dsint(-ha,ra,0.,0.,0.,0.,ra) + eta*dsint(-ha,ra,2*ha,0.,0.,0.,ra) - (1+eta)*dsint(-ha,ra,ha+hb,0.,-dy,0.,ra) + eta*dsint(-ha,ra,0.,2*ha,0.,0.,ra) + eta*eta*dsint(-ha,ra,2*ha,2*ha,0.,0.,ra) - eta*(1+eta)*dsint(-ha,ra,ha+hb,2*ha,-dy,0.,ra) - (1+eta)*dsint(-ha,ra,0.,ha+hb,0.,-dy,ra) - eta*(1+eta)*dsint(-ha,ra,2*ha,ha+hb,0.,-dy,ra) + (1+eta)*(1+eta)*dsint(-ha,ra,ha+hb,ha+hb,-dy,-dy,ra);
valA1 /= 8*PI*e1;
valA2 = dsint(-ra,-ha,ha+hb,ha+hb,-dy,-dy,ra) - eta*dsint(-ra,-ha,ha-hb,ha+hb,-dy,-dy,ra) - (1-eta)*dsint(-ra,-ha,0.,ha+hb,0.,-dy,ra) - eta*dsint(-ra,-ha,ha+hb,ha-hb,-dy,-dy,ra) + eta*eta*dsint(-ra,-ha,ha-hb,ha-hb,-dy,-dy,ra) + eta*(1-eta)*dsint(-ra,-ha,0.,ha-hb,0.,-dy,ra) - (1-eta)*dsint(-ra,-ha,ha+hb,0.,-dy,0.,ra) + eta*(1-eta)*dsint(-ra,-ha,ha-hb,0.,-dy,0.,ra) + (1-eta)*(1-eta)*dsint(-ra,-ha,0.,0.,0.,0.,ra);
valA2 /= 8*PI*e2;
valB1 = dsint(hb,rb,-ha-hb,-ha-hb,dy,dy,rb) + eta*dsint(hb,rb,ha-hb,-ha-hb,dy,dy,rb) - (1+eta)*dsint(hb,rb,0.,-ha-hb,0.,dy,rb) + eta*dsint(hb,rb,-ha-hb,ha-hb,dy,dy,rb) + eta*eta*dsint(hb,rb,ha-hb,ha-hb,dy,dy,rb) - eta*(1+eta)*dsint(hb,rb,0.,ha-hb,0.,dy,rb) - (1+eta)*dsint(hb,rb,-ha-hb,0.,dy,0.,rb) - eta*(1+eta)*dsint(hb,rb,ha-hb,0.,dy,0.,rb) + (1+eta)*(1+eta)*dsint(hb,rb,0.,0.,0.,0.,rb);
valB1 /= 8*PI*e1;
valB2 = dsint(-rb,hb,0.,0.,0.,0.,rb) - eta*dsint(-rb,hb,-2*hb,0.,0.,0.,rb) - (1-eta)*dsint(-rb,hb,-ha-hb,0.,dy,0.,rb) - eta*dsint(-rb,hb,0.,-2*hb,0.,0.,rb) + eta*eta*dsint(-rb,hb,-2*hb,-2*hb,0.,0.,rb) + eta*(1-eta)*dsint(-rb,hb,-ha-hb,-2*hb,dy,0.,rb) - (1-eta)*dsint(-rb,hb,0.,-ha-hb,0.,dy,rb) + eta*(1-eta)*dsint(-rb,hb,-2*hb,-ha-hb,0.,dy,rb) + (1-eta)*(1-eta)*dsint(-rb,hb,-ha-hb,-ha-hb,dy,dy,rb);
valB2 /= 8*PI*e2;
val = valA1 + valA2 + valB1 + valB2;
break;
}
return(val);
}
double dsint(zmin,zmax,alphaz,betaz,alphay,betay,rr) double zmin,zmax,alphaz,betaz,alphay,betay,rr;
{
double lower, upper, value, gauss(), func();
ya = alphay;
yb = betay;
za = alphaz;
zb = betaz;
r = rr;
lower = zmin;
upper = zmax;
value = gauss(NSUB,lower,upper,func);
return(value);
}
double func(z) double z;
{
double lower1, upper1, f, gauss1(), func1();
zf = z;
lower1 = 0.;
upper1 = twoPI; f = gauss1(NSUB,lower1,upper1,func1);
return (f);
}
double func1(y) double y;
{
double f1,d1,d2;
f1 = r*r+sqrt(r*r-zf*zf)*sin(y)*yb+zf*zb;
d1 = r*r+ya*ya+za*za+2*ya*sqrt(r*r-zf*zf)*sin(y)+2*za*zf;
d2 = r*r+yb*yb+zb*zb+2*yb*sqrt(r*r-zf*zf)*sin(y)+2*zb*zf;
f1 /= (sqrt(d1*d2)*d2);
return (f1);
}
double gauss(nsub,a,b,func) double a, b, (*func)();
int nsub;
{
static double legend[] = {
0.095012509837637, 0.281603550779258, 0.458016777657227, 0.617876244402643, 0.755404408355003, 0.865631202387831, 0.944575023073232, 0.989400934991649}
;
static double weight[] = {
0.189450610455068, 0.182603415044923, 0.169156519395002, 0.149595988816576, 0.124628971255533, 0.095158511682492, 0.062253523938647, 0.027152459411754}
;
double next_a, dx, fy, yplus, yminus, integral = 0.;
int i, j;
dx = (b-a)/nsub; /* size of each subinterval */
next_a = a;
for (j=0; j<= nsub-1; j++) {
a = next_a;
b = a + dx;
next_a = b;
for (i=0; i<= 7; i++) /* The 16-point Gaussian quadrature loop */ {
yplus = ( dx * legend[i] + b + a)/2.;
yminus = ( -dx * legend[i] + b + a)/2.;
fy = (*func) (yplus) + (*func) (yminus);
integral += weight[i]*fy;
}
}
/* end loop over the nsub sub_intervals */ integral *= dx/2.;
return(integral);
}
double gauss1(nsub,a,b,func1) double a, b, (*func1)();
int nsub;
{
static double legend1[] = {
0.095012509837637, 0.281603550779258, 0.458016777657227, 0.617876244402643, 0.755404408355003, 0.865631202387831, 0.944575023073232, 0.989400934991649}
;
static double weight1[] = {
0.189450610455068, 0.182603415044923, 0.169156519395002, 0.149595988816576, 0.124628971255533, 0.095158511682492, 0.062253523938647, 0.027152459411754}
;
double next_a, dx, fy, yplus, yminus, integral = 0.;
int i, j;
dx = (b-a)/nsub; /* size of each subinterval */
next_a = a;
for (j=0; j<= nsub-1; j++) {
a = next_a;
b = a + dx;
next_a = b;
for (i=0; i<= 7; i++) /* The 16-point Gaussian quadrature loop */ {
yplus = ( dx * legend1[i] + b + a)/2.;
yminus = ( -dx * legend1[i] + b + a)/2.;
fy = (*func1) (yplus) + (*func1) (yminus);
integral += weight1[i]*fy;
}
}
/* end loop over the nsub sub_intervals */ integral *= dx/2.;
return(integral);
}