subroutine trin(x1,y1,z1,x2,y2,z2,x3,y3,z3,a,pp,xp1,yp1,zp1, 1 xp2,yp2,zp2,xp3,yp3,zp3,c1,c2,c3) c c********************************************************************** c (1) function : c generate the transformation matrix between global coord. and c local coord. c********************************************************************** c real*8 a(3,3),xg(4),yg(4),xpg(4),ypg(4),zpg(4) real*8 x1,y1,z1,x2,y2,z2,x3,y3,z3,pp,xp1,yp1,zp1,xp2,yp2,zp2 real*8 xp3,yp3,zp3,c1,xl1,xl2,xl3,c2,c3,b1,b2,b3 real*8 a1,a2,a3,p1,p2,p3 real*8 ksi,eta,zta,f common/ggp/ xpg,ypg,zpg f(x1,y1,z1,x2,y2,z2)=dsqrt((x1-x2)**2+(y1-y2)**2+(z1-z2)**2) xl1=f(x2,y2,z2,x3,y3,z3) xl2=f(x1,y1,z1,x3,y3,z3) xl3=f(x1,y1,z1,x2,y2,z2) c** c** (a1,a2,a3) : unit vector in x'-direction c** a1=(x1-x3)/xl2 a2=(y1-y3)/xl2 a3=(z1-z3)/xl2 c** c** (c1,c2,c3) : unit vector in z'-direction (normal) c** p1=(y1-y3)*(z2-z3)-(y2-y3)*(z1-z3) p2=(x2-x3)*(z1-z3)-(x1-x3)*(z2-z3) p3=(x1-x3)*(y2-y3)-(x2-x3)*(y1-y3) pp=dsqrt(p1*p1+p2*p2+p3*p3) c1=p1/pp c2=p2/pp c3=p3/pp c** c** (b1,b2,b3) : unit vector in y'-direction c** b1=c2*a3-c3*a2 b2=c3*a1-c1*a3 b3=c1*a2-c2*a1 c** c** a(i,j) : transformation matrix c** a(1,1)=b2*c3-c2*b3 a(1,2)=c1*b3-b1*c3 a(1,3)=b1*c2-c1*b2 a(2,1)=c2*a3-a2*c3 a(2,2)=a1*c3-c1*a3 a(2,3)=c1*a2-a1*c2 a(3,1)=a2*b3-b2*a3 a(3,2)=b1*a3-a1*b3 a(3,3)=a1*b2-b1*a2 c** c** (xp1,yp1,zp1) : new coordinates of point(x1,y1,z1) c** ksi=x1-x3 eta=y1-y3 zta=z1-z3 call new(a,ksi,eta,zta,xp1,yp1,zp1) c** c** (xp2,yp2,zp2) : new coordinates of point(x2,y2,z2) ksi=x2-x3 eta=y2-y3 zta=z2-z3 call new(a,ksi,eta,zta,xp2,yp2,zp2) c** c** (xp3,yp3,zp3) : new coordinates of point(x3,y3,z3) c** xp3=0.0 yp3=0.0 zp3=0.0 call gaussp(xp1,xp2,yp2,pp,xg,yg) do 10 i=1,3 xpg(i)=x3+a1*xg(i)+b1*yg(i) ypg(i)=y3+a2*xg(i)+b2*yg(i) zpg(i)=z3+a3*xg(i)+b3*yg(i) 10 continue c write(*,*) 'x1,y1,z1',x1,y1,z1 c write(*,*) 'x2,y2,z2',x2,y2,z2 c write(*,*) 'x3,y3,z3',x3,y3,z3 c do 20 i=1,4 c ksi=xpg(i)-x3 c eta=ypg(i)-y3 c zta=zpg(i)-z3 c call new(a,ksi,eta,zta,xx,yy,zz) c write(*,*) i,xpg(i),ypg(i),zpg(i),xx,yy,zz c20 continue return end subroutine gaussp(xp1,xp2,yp2,pp,xg,yg) real*8 xg(4),yg(4),xp1,xp2,yp2,pp real*8 l1(3),l2(3),l3(3) c data l1/0.333333,0.6,0.2,0.2/,l2/0.333333,0.2,0.6,0.2/ data l1/0.5,0.0,0.5/,l2/0.5,0.5,0.0/ do 10 i=1,3 xg(i)=pp*(l1(i)+xp2*l2(i)/xp1)/yp2 yg(i)=pp*l2(i)/xp1 10 continue c write(*,*) 'xp1,xp2,yp2,pp',xp1,xp2,yp2,pp do 20 i=1,3 l3(i)=(xp1*yp2-yp2*xg(i)+(xp2-xp1)*yg(i))/pp c write(*,*) 'l1,l2,l3',l1(i),l2(i),l3(i) c write(*,*) i,xg(i),yg(i) 20 continue return end subroutine corr(xt1,yt1,zt1,xt2,yt2,zt2,xt3,yt3,zt3,l, 1 x1,y1,z1,x2,y2,z2,x3,y3,z3) real*8 xt1,yt1,zt1,xt2,yt2,zt2,xt3,yt3,zt3,x1,y1,z1,x2,y2,z2 real*8 x3,y3,z3 if(l.eq.1) go to 10 if(l.eq.2) go to 40 if(l.eq.3) go to 30 if(l.eq.4) go to 20 10 x1=xt1 y1=yt1 z1=zt1 x2=xt2 y2=yt2 z2=zt2 x3=xt3 y3=yt3 z3=zt3 go to 50 20 x1=xt1 y1=-yt1 z1=zt1 x2=xt3 y2=-yt3 z2=zt3 x3=xt2 y3=-yt2 z3=zt2 go to 50 30 x1=xt1 y1=-yt1 z1=-zt1 x2=xt2 y2=-yt2 z2=-zt2 x3=xt3 y3=-yt3 z3=-zt3 go to 50 40 x1=xt1 y1=yt1 z1=-zt1 x2=xt3 y2=yt3 z2=-zt3 x3=xt2 y3=yt2 z3=-zt2 50 return end subroutine cacu(xc,yc,zc,x1,y1,z1,x2,y2,z2,x3,y3,z3,h,g) c c************************************************************************ c (1) function : c compute potential at (xc,yc,zc),induced by unit source (doublet) c triangle (x1,y1,z1),(x2,y2,z2) and (x3,y3,z3) c (2) variable : c h : potntial induced by unit doublet triangle c g : potntial induced by unit source triangle c************************************************************************ c real*8 xc,yc,zc,x1,y1,z1,x2,y2,z2,x3,y3,z3,h,g real*8 xp1,yp1,xp2,yp2,xp3,yp3,xs1,xs2,tempx,tempy,h1,g1 real*8 h2,g2 1 xp1=x1-xc yp1=y1-yc xp2=x2-xc yp2=y2-yc xp3=x3-xc yp3=y3-yc xs1=xp1-yp1*(xp2-xp1)/(yp2-yp1) xs2=xp2-yp2*(xp2-xp3)/(yp2-yp3) c if(nk.eq.398.and.nh.eq.445) print *,"xs1=",xs1," xs2=",xs2 if(yp1*yp2.lt.0..and.xs1.lt.0..and.xs2.lt.0.) go to 10 go to 20 10 xp1=-xp1 xp2=-xp2 xp3=-xp3 tempx=xp1 tempy=yp1 xp1=xp3 yp1=yp3 xp3=tempx yp3=tempy c if(nk.eq.398.and.nh.eq.445) print *,xp1,yp1,xp2,yp2,xp3,yp3 20 xp=0. yp=0. call intt(xp1,yp1,xp2,yp2,zc,h1,g1) call intt(xp2,yp2,xp3,yp3,zc,h2,g2) h=h1+h2 g=g1+g2 return end subroutine new(a,ksi,eta,zta,xp,yp,zp) c c************************************************************************* c (1) function : c transform global coord. (ksi,eta,zta) to local coord. (xp,yp,zp) c (2) variables : c a : transformation matrix c************************************************************************* c real*8 a(3,3) real*8 ksi,eta,zta,xp,yp,zp xp=ksi*a(1,1)+eta*a(1,2)+zta*a(1,3) yp=ksi*a(2,1)+eta*a(2,2)+zta*a(2,3) zp=ksi*a(3,1)+eta*a(3,2)+zta*a(3,3) return end subroutine intt(x1,y1,x2,y2,zp,h,g) c c************************************************************************** c (1) function : c line integral along edges of the triangle c (2) criteria for convergence of intgegral : 1.e-2 c (h2-h1)/h2 < 1.e-2 c (g2-g1)/g2 < 1.e-2 c where c h1,g1 : numerical value of first integral c h2,g2 : numerical value of second integral c************************************************************************** c real*8 x1,y1,x2,y2,zp,h,g,delx,dely,ddx,ddy,h2,g2,xx1,yy1,xx2,yy2 real*8 hh,gg,ch1,cg1 delx=x2-x1 dely=y2-y1 do 10 i=1,20 ddx=delx/i ddy=dely/i h2=0.0 g2=0.0 do 20 j=1,i xx1=x1+(j-1)*ddx yy1=y1+(j-1)*ddy xx2=x1+j*ddx yy2=y1+j*ddy call gauss(xx1,yy1,xx2,yy2,zp,hh,gg) h2=h2+hh g2=g2+gg 20 continue if(i.eq.1) go to 30 if(abs(h2).lt.1.e-8) go to 50 ch1=(h2-h1)/h2 cg1=(g2-g1)/g2 if(abs(ch1).le.1.e-5.and.abs(cg1).le.1.e-5) go to 40 go to 30 50 cg1=(g2-g1)/g2 if(abs(cg1).le.1.e-5) go to 40 30 h1=h2 g1=g2 10 continue 40 h=h2 g=g2 return end subroutine gauss(x1,y1,x2,y2,zp,h,g) c c************************************************************************** c gaussin quadrature (five points) c variables : c gn : roots c w : weight for each root c************************************************************************** c real*8 gn(5),w(5),xx(5),yy(5),q real*8 delx,dely,s,cs,x1,y1,x2,y2,h,g,dd,xl real*8 xq,yq,qq,zp,d data gn/-0.9061798459,-0.5384693101,0.0, 1 0.5384693101, 0.9061798459/ data w/ 0.2369268851, 0.4786286705, 0.5688888889, 1 0.4786286705, 0.2369268851/ delx=x2-x1 dely=y2-y1 if(abs(delx).le.1.e-10) go to 100 s=dely/delx cs=dsqrt(s*s/(1.+s*s)) if(y2.lt.y1) cs=-cs c print *,"cs=",cs xl=abs(delx)*sqrt(1.+s*s) dd=2.*delx*(1.+s*s) do 10 i=1,5 xx(i)=xl*xl*(gn(i)+1)/dd+x1 yy(i)=(xx(i)-x1)*s+y1 10 continue go to 40 100 xl=abs(dely) cs=1. if(y2.lt.y1) cs=-cs do 110 i=1,5 xx(i)=x1 yy(i)=(gn(i)+1.)/2.*(y2-y1)+y1 110 continue 40 h=0. g=0. do 120 i=1,5 xq=xx(i) yq=yy(i) qq=q(xq,yq,zp,d) c print*,"xq=",xq," yq=",yq," hh=",qq," gg=",d h=h+qq*w(i) g=g+d*w(i) 120 continue h=h*cs*xl/2. g=g*cs*xl/2. c print *,"h=",h," g=",g return end function q(xq,yq,zp,d) real*8 q,xq,yq,zp,d,a,b,dd a=yq*yq+zp*zp b=dsqrt(xq*xq+yq*yq+zp*zp) if(zp.eq.0.) q=0. if(zp.ne.0.) q=-(zp*xq)/(a*b) dd=xq+b d=-dlog(xq+b) return end