san.zip.txt

求三极法井间低阻球体上视电阻率值的程序

/*跨孔条件下，井中10m供电二极剖面法球体程序（球体中心（坐标原点）在两井中心的左侧）*/
#include <stdio.h>
#include <math.h>
double P(int N, double X)
{   

   if(N == 0)
       return 1;
   else
       if(N== 1)  
           return X;
       else
           return ((2*N-1)*X*P(N-1, X) - (N-1) * P(N-2, X))/N;                                    
}                                                 

double Sum(double r0 ,double p1,double p2,double h0,double x,double hA,double hM)

{ int n;  double a, b, c,d1;double sum=0.0; double f,k,Ha,Hm,C1,R=78,D1;/* d1=h0/r0，C1为d1/r0，D1为r/r0 ,Ha=hA/r0,Hm=hM/r0*/
for(n=0;n<8;n++)
	{ a=R/r0;b=p2/p1;c=x/r0;d1=h0/r0;Ha=hA/r0;Hm=hM/r0;
  C1=sqrt((0.5*a-c)*(0.5*a-c)+(d1-Ha)*(d1-Ha));D1=sqrt((0.5*a+c)*(0.5*a+c)+(d1-Hm)*(d1-Hm));
       k=(D1*D1+C1*C1-a*a)/(2*D1*C1);
       f=pow(D1,-(n+1))*pow(C1,-(n+1));
       sum=sum+n*(b-1)/(n+(n+1)*b)*f*P(n, k);
	}

 return sum ; 
}


double y()
{double l=78,r0=20 ,p2=50,p1=100,h0=20,x=1,hA=10,I=10 ,y1,U,AM;/* a=R/ro，b=P2/P1，c=x/r0，d1=h0/r0；，其中，R为AM之间的距离，r0为球体半径，h0为球体埋深，P1，P2分别为围岩与球体的电阻率值。x为oo’的距离，即球心相对两井中点的水平距离。o为求球心在地面上的投影点，o'为AB两钻孔的中点*/
static double hM=0;
AM=sqrt(l*l+(hM-hA)*(hM-hA));
y1= 1+AM/r0*Sum(r0 ,p1,p2,h0,x,hA,hM);/*y为二极装置下的理论视电阻率值*/
 
U=I/(4*3.1416*AM)*y1;
hM=hM+0.5;
 return y1;
 	}
void main()

{ double l=78,r0=20 ,p2=50,p1=100,h0=20,x=1,hA=10,MN=1,hn,I=10 ,y2,y3,U,R,AM;/* a=R/ro，b=P2/P1，c=x/r0，d1=h0/r0；，其中，R为AM之间的距离，r0为球体半径，h0为球体埋深，P1，P2分别为围岩与球体的电阻率值。x为oo’的距离，即球心相对两井中点的水平距离。o为求球心在地面上的投影点，o'为AB两钻孔的中点*/

  for(hn=1;hn<=40.5;hn=hn+0.5) 

{ R=sqrt(l*l+(hn-hA)*(hn-hA));
 
y2= 1+R/r0*Sum(r0 ,p1,p2,h0,x,hA,hn);/*y为二极装置下的理论视电阻率值*/
 U=I/(4*3.1416*R)*y2;
static double hM=0;
AM=sqrt(l*l+(hM-hA)*(hM-hA));
 y3= 1.0/MN*(R*y()-AM*y2);   

 printf("%f %f\n", y3,U);

hM=hM+0.5;
 } 

}