［Okada(1985)の手法検証プログラム］

      dimension disp(3)
c
      rmu=1.0
      ramda=1.0
      xx=2.0
      yy=3.0
      dd=4.0
      dl=3.0
      dw=2.0
      u1=0.0
      u2=0.0
      u3=1.0
      delta=70.0
c
      call okada85(disp,rmu,ramda,xx,yy,dd,delta,dl,dw,u1,u2,u3)
c
      write(*,*) 'Ux=',disp(1)
      write(*,*) 'Uy=',disp(2)
      write(*,*) 'Uz=',disp(3)
c
      stop
      end
c **********************************************************************
      subroutine okada85(disp,rmu,ramda,xx,yy,dd,delta,dl,dw,u1,u2,u3)
c **********************************************************************
      dimension disp(3),o1(3),o2(3),o3(3),o4(3)
c
      pai=2.0*acos(0.0)
      p=yy*cos(pai*delta/180.0)+dd*sin(pai*delta/180.0)
c
      gzi=xx
      rita=p
      call okadaff(rmu,ramda,yy,gzi,rita,dd,delta,u1,u2,u3,o1)
      gzi=xx
      rita=p-dw
      call okadaff(rmu,ramda,yy,gzi,rita,dd,delta,u1,u2,u3,o2)
      gzi=xx-dl
      rita=p
      call okadaff(rmu,ramda,yy,gzi,rita,dd,delta,u1,u2,u3,o3)
      gzi=xx-dl
      rita=p-dw
      call okadaff(rmu,ramda,yy,gzi,rita,dd,delta,u1,u2,u3,o4)
c
      disp(1)=o1(1)-o2(1)-o3(1)+o4(1)
      disp(2)=o1(2)-o2(2)-o3(2)+o4(2)
      disp(3)=o1(3)-o2(3)-o3(3)+o4(3)
c
      return
      end
c **********************************************************************
      subroutine okadaff(rmu,ramda,yy,gzi,rita,dd,delta,u1,u2,u3,oo)
c **********************************************************************
      dimension oo(3)
c
      pai=2.0*acos(0.0)
      delsin=sin(pai*delta/180.0)
      delcos=cos(pai*delta/180.0)
      p=yy*delcos+dd*delsin
      q=yy*delsin-dd*delcos
      y2=rita*delcos+q*delsin
      d2=rita*delsin-q*delcos
      rr=sqrt(gzi**2.0+rita**2.0+q**2.0)
      xx=sqrt(gzi**2.0+q**2.0)
      rmm=rmu/(ramda+rmu)
c
      if(abs(delta).eq.90.0 .or. abs(delta).eq.270.0) then
        delcos=0.0
        ri5=-1.0*rmm*gzi*delsin/(rr+d2)
        ri4=-1.0*rmm*q/(rr+d2)
        ri3=0.5*rmm*(rita/(rr+d2)+y2*q/(rr+d2)**2.0-alog(rr+rita))
        ri2=-1.0*rmm*alog(rr+rita)-ri3
        ri1=-0.5*rmm*gzi*q/(rr+d2)**2.0
      else
        ri5=rmm*2.0/delcos*atan((rita*(xx+q*delcos)+xx*(rr+xx)*delsin)
     &           /(gzi*(rr+xx)*delcos))
        ri4=rmm/delcos*(alog(rr+d2)-delsin*alog(rr+rita))
        ri3=rmm*(y2/(delcos*(rr+d2))-alog(rr+rita))+delsin*ri4/delcos
        ri2=-1.0*rmm*alog(rr+rita)-ri3
        ri1=rmm*(-1.0*gzi/(delcos*(rr+d2)))-delsin*ri5/delcos
      end if
c
      ux1=-0.5*u1/pai*(gzi*q/(rr*(rr+rita))+atan(gzi*rita/(q*rr))
     &                                                +ri1*delsin)
      uy1=-0.5*u1/pai*(y2*q/(rr*(rr+rita))+q*delcos/(rr+rita)
     &                                                +ri2*delsin)
      uz1=-0.5*u1/pai*(d2*q/(rr*(rr+rita))+q*delsin/(rr+rita)
     &                                                +ri4*delsin)
      ux2=-0.5*u2/pai*(q/rr-ri3*delsin*delcos)
      uy2=-0.5*u2/pai*(y2*q/(rr*(rr+gzi))
     &           +delcos*atan(gzi*rita/(q*rr))-ri1*delsin*delcos)
      uz2=-0.5*u2/pai*(d2*q/(rr*(rr+gzi))
     &           +delsin*atan(gzi*rita/(q*rr))-ri5*delsin*delcos)
      ux3=0.5*u3/pai*(q**2.0/(rr*(rr+rita))-ri3*delsin**2.0)
      uy3=0.5*u3/pai*(-1.0*d2*q/(rr*(rr+gzi))
     &          -delsin*(gzi*q/(rr*(rr+rita))-atan(gzi*rita/(q*rr)))
     &                                              -ri1*delsin**2.0)
      uz3=0.5*u3/pai*(y2*q/(rr*(rr+gzi))
     &          +delcos*(gzi*q/(rr*(rr+rita))-atan(gzi*rita/(q*rr)))
     &                                              -ri5*delsin**2.0)
c
      oo(1)=ux1+ux2+ux3
      oo(2)=uy1+uy2+uy3
      oo(3)=uz1+uz2+uz3
c
      return
      end