program v2fid c*********************************************************************** c Deformation due to Dislocation Sources on a Finite Fault c CONDITIONS: c Structure: elastic-Viscoelastic, 2 Layer c Source: Finite fault, arbitrary slip c Deformation: Internal displacement, all components c for various Times c Coordinate system is the same as Fukahata & Matsu'ura (2005,GJI) c The strike of finite fault is taken to be the x-axis. c (i.e. Observation points align parallel to the fault strike) c c modified from v2faiat_Isa2.f and v2faiatAo.f c coded at 050815, modified at 060630, 151120, 170829 c IMPORTANT!! Correction to the output comment (170829): c "ux, uy, uz [mm]" ---> "uz, ux, uy [mm]" c c References: Fukahata & Matsu'ura, 2005, GJI, 161, 507-521. c Fukahata & Matsu'ura, 2006, GJI, 166, 418-434. c c## The program is coded by Yukitoshi Fukahata. c*********************************************************************** implicit double precision(a-h,o-z) c Paramters that begin with a-h and o-z are double precision, c and that begin with i-n are integer. parameter (nt=2, pi=3.14159265358979, eps=1.0D-7) character outputfile*20 CHANGE1/3 ************************************************************** c ### Parameters for accuracy in computation ### parameter (pr_gze=18.0, gzs=0.002, zdmin=0.05) c ### Output File Name ### parameter (outputfile='temp.dat') c ### Source Parameters ### parameter (du=1000, id1=10,id2=15, L=10, ded=45.0, chid=60.0) parameter (isx=1, isz=1) c ### Parameters for Observation ### parameter (nti=2) parameter (nxi=-5,nxe=15,nxs=10, nyi=-3,nye=9,nys=6, & nzi=0,nzn=2,nzs=40) c NECESSARY CONDITIONS FOR PARAMETERS -------------------- c ded =! 0 c (id2-id1)/isz = integer, L/isx = integer, nxs/isx = integer c (L+nxe-nxi)*(nye-nyi)*nzn*nti should not be very big. c (Minimum distance from fault to ob. point) > 4*isx is preferable c*********************************************************************** parameter (npxi=nxi-L) dimension sux(nxi:nxe,nyi:nye,nzn,nti), & suy(nxi:nxe,nyi:nye,nzn,nti), & suz(nxi:nxe,nyi:nye,nzn,nti) dimension ux(npxi:nxe,nyi:nye,nzn,nti), & uy(npxi:nxe,nyi:nye,nzn,nti), & uz(npxi:nxe,nyi:nye,nzn,nti) c----------------------------------------------------------------------- c unit of length: [km] c nt : number of layers c pr_gze: cut off wavenumber c gzs : step of integral variable (wavenumber) "gz" c zdmin : minimum depth difference between source and observation c outputfile: subsutitute the name of output file c c -------------- Source parameters -------------------- c du : dislocation on the fault [mm] c id1 : upper bound of the fault [km] c id2 : lower bound of the fault [km] c L : horizontal length of the fault [km] c ded : delta, dip of the fault [degree] c chid: slip angle (rake) of the fault [degree] c Fault Geometry (coordinates of apexes): [(0,0,id1), (L,0,id1), c (L,(id2-id1)/tan(ded),id2), (0,(id2-id1)/tan(ded),id2)] c isx (isz): horizontal (vertical) span of each fault patch [km] c c --- in the computation --- c d : depth of a source [km] c de : delta, dip of the fault [radian] c chi: slip angle (rake) of the fault [radian] c fai: fai, angle from the strike of the fault (see Fig.1b) [radian] c ys : y-coordinate of a source [km] c dd : area of each fault patch [km^2] c ----------------------------------------------------- c c ---------- Parameters for observation points -------- c nti: number of times to calculate deformation c nxi -- nxe: calculation range of x [km] c nyi -- nye: calculation range of y [km] c nxs (nys) : grid interval in x (y) direction [km] c for "z" : The program calculate nzn times for z c from nzi with a regular interval nzs [km] as below c do 55 nz=1,nzn c z = dble(nzi+(nz-1)*nzs) c c ux,uy,uz : displacement due to a point dislocation on each patch c =G(x,y,z; xs,ys,d) [mm] c sux,suy,suz: displacement after taking summation on x [mm] c----------------------------------------------------------------------- dimension Ytot(4,3,nti),Ydtot(2,2,nti) dimension af(3,3),afd(2,2), besr(3,3),besrd(2,2) dimension vp(nt),vs(nt), ti(nti) common /B1/ gam(2) common /B2/ um(2) common /B3/ rou(2) common /B4/ gra common /C1/ HL common /D1/ ak common /D2/ tau CHANGE2/3 ********* Times and Structural Parameters ******************** data ti/0.d0, 10.d0/ data vp/6.0d0, 8.0d0/ vs/3.5d0, 4.5d0/ data rou/2.6d12, 3.4d12/ data gra/9.8d-3/ c data gra/0.0d0/ data vis/1.0d19/ data HL/30.d0/ c*********************************************************************** c---------------------------------------------------------------- c ti : times to calculate deformation [yr or tau] c (tau: Maxwell relaxation time) c vp : P-wave velocity [km/s] c vs : S-wave velocity [km/s] c gam: gamma c um : rigidity [Pa']=[N'/m^2']=[kg km/s^2]*[1/km^2]=[kg/(km s^2)] c rou: density [kg/km^3] c gra: gravity accerelation [km/s^2] c vis: viscosity [Pa s] c HL : depth of the layer interface [km] c---------------------------------------------------------------- c --------------- Initialization ------------------- open(10,file=outputfile) do 10 i=1,nt gam(i) = 1-(vs(i)/vp(i))**2 um(i) = rou(i)*vs(i)*vs(i) 10 continue visd = vis*1000. CHANGE 3/3 ************ Scale of time (year or tau) ******************** cc time in year tau = visd/um(2)/3.1556926d7 cc time in tau c tau=1. c*********************************************************************** write(*,*) 'tau=', real(tau) ak = (4*gam(2)-1)/3. de = ded/180.0*pi chi = chid/180.0*pi dd = dble(isx*isz)/sin(de) if(abs(de).lt.eps) write(*,*) 'error de' if(mod(nxs,isx).ne.0) write(*,*) 'nxs/isx is not integer' do 19 it=1,nti do 20 nz=1,nzn do 21 nx=npxi,nxe,isx do 22 ny=nyi,nye,nys uz(nx,ny,nz,it)=0. ux(nx,ny,nz,it)=0. uy(nx,ny,nz,it)=0. 22 continue 21 continue 20 continue 19 continue c ------------------ Main Calculation ------------------ w1 = gzs/4./pi*du*dd njz = (id2-id1)/isz do 50 jz = 1,njz wz=(jz-0.5)*isz d = id1+wz call wh_Lay(d,m) if(cos(de).lt.eps) then ys=0. else ys=wz/tan(de) endif write(*,*) jz, real(ys),real(d),real(dd),real(de),m do 55 nz=1,nzn z = dble(nzi+(nz-1)*nzs) zd = z-d if(zd.le.zdmin .and. zd.ge.0.) then c if(nz.eq.int(Hl)) then c z=d-zdmin c else z=d+zdmin c endif write(*,*) 'i,nz,z,d',i,nz,z,d,' +' else if(zd.ge.-zdmin .and. zd.lt.0.) then if(d.lt.zdmin) then z=d+zdmin else z=d-zdmin endif write(*,*) 'i,nz,z,d',i,nz,z,d,' -' end if call wh_Lay(z,j) c########## For check of the kernel program ################# c gz = 0.08 c write(*,*) 'z,',z,j,m,d c call mYdtot(z,j,m,d,gz, nti,ti, Ydtot) c call mYtot(z,j,m,d,gz, nti,ti, Ytot) c write(*,*) 'z,U1',real(z),j,(real(Ytot(1,k,1)),k=1,3) c write(*,*) 'z,U2',real(z),j,(real(Ytot(2,k,1)),k=1,3) c write(*,*) (real(Ydtot(1,k,1)),k=1,2) c 55 continue c 50 continue c stop c############################################################ gze = pr_gze/abs(z-d) ngze= int(gze/gzs) do 100 igz=1,ngze gz = igz*gzs - gzs/2. call mYtot(z,j,m,d,gz, nti,ti, Ytot) call mYdtot(z,j,m,d,gz, nti,ti, Ydtot) do 110 ny=nyi,nye,nys y = ny - ys do 115 nx=npxi+isx,nxe,isx x = nx-0.5*isx r = dsqrt(y*y + x*x) if((abs(y).lt.eps) .and. (abs(x).eq.eps)) then fai=pi/2. else fai = datan2(y,x) endif call vecaf(de,fai,chi, af,afd) call vbess(r,gz, besr,besrd) do 119 it=1,nti wur=0. wuf=0. do 120 k=1,3 wur = wur + Ytot(1,k,it)*af(1,k)*besr(1,k) wuf = wuf + Ytot(1,k,it)*af(2,k)*besr(2,k) uz(nx,ny,nz,it) = uz(nx,ny,nz,it) & + Ytot(2,k,it)*af(3,k)*besr(3,k) 120 continue do 125 k=1,2 wur = wur + Ydtot(1,k,it)*afd(1,k)*besrd(1,k) wuf = wuf - Ydtot(1,k,it)*afd(2,k)*besrd(2,k) 125 continue ux(nx,ny,nz,it) = ux(nx,ny,nz,it) & + (cos(fai)*wur - sin(fai)*wuf) uy(nx,ny,nz,it) = uy(nx,ny,nz,it) & + (sin(fai)*wur + cos(fai)*wuf) 119 continue 115 continue 110 continue 100 continue 55 continue 50 continue c ----------------- Output ---------------- write(10,*) 'time[yr or tau]; z, x, y [km]; uz, ux, uy [mm] (upl &ift positive in uz)' do 200 it=1,nti do 200 nz=1,nzn do 210 ny=nyi,nye,nys do 220 nx=nxi,nxe,nxs sux(nx,ny,nz,it)=0. suy(nx,ny,nz,it)=0. suz(nx,ny,nz,it)=0. do 230 nsx=nx-L+isx,nx,isx sux(nx,ny,nz,it)=sux(nx,ny,nz,it) + ux(nsx,ny,nz,it)*w1 suy(nx,ny,nz,it)=suy(nx,ny,nz,it) + uy(nsx,ny,nz,it)*w1 suz(nx,ny,nz,it)=suz(nx,ny,nz,it) + uz(nsx,ny,nz,it)*w1 230 continue write(10,1010) real(ti(it)),nzi+(nz-1)*nzs, nx,ny, $ -real(suz(nx,ny,nz,it)), $ real(sux(nx,ny,nz,it)),real(suy(nx,ny,nz,it)) 1010 format(F10.3, 3I6, 3E14.5) 220 continue 210 continue 200 continue end c/////////////////////////////////////////////////////////////////////// subroutine wh_Lay(x,n) c-------------------------------------------------------- c which(n-th) Layer does x belong to ? c-------------------------------------------------------- implicit double precision(a-h,o-z) common /C1/ HL if ((x.gt.-0.0001).and.(x.le.HL)) then n=1 else if (x.gt.HL) then n=2 else write(*,*) 'wh_Lay is error!' stop endif end c----------------------------------------------------------------------- subroutine vecaf(de,fai,chi, af,afd) c-------------------------------------------------------- c vecaf: vector a(fai) c-------------------------------------------------------- implicit double precision(a-h,o-z) dimension af(3,3),afd(2,2) sx=sin(chi) cx=cos(chi) d2=2*de f2=2*fai af(1,1) = sx*sin(d2)/4. af(1,2) =-sx*cos(d2)*sin(fai) + cx*cos(de)*cos(fai) af(1,3) = sx*sin(d2)*cos(f2)/4. + cx*sin(de)*sin(f2)/2. af(2,1) = 0. af(2,2) =-sx*cos(d2)*cos(fai) - cx*cos(de)*sin(fai) af(2,3) =-sx*sin(d2)*sin(f2)/2.+ cx*sin(de)*cos(f2) af(3,1) = af(1,1) af(3,2) = af(1,2) af(3,3) = af(1,3) afd(1,1) = -af(1,2) afd(1,2) = -4*af(1,3) afd(2,1) = af(2,2) afd(2,2) = af(2,3) end c----------------------------------------------------------------------- subroutine vbess(r,gz, besr,besrd) c--------------------------------------------------------- c vbess: vector J(r,gz) c--------------------------------------------------------- implicit double precision(a-h,o-z) dimension bes(0:3), besr(3,3),besrd(2,2) x=r*gz bes(0) = bessj0(x) bes(1) = bessj1(x) bes(2) = bessj(2,x) bes(3) = bessj(3,x) besr(1,1) =-gz*bes(1) besr(1,2) = gz*(bes(0)-bes(2))/2. besr(1,3) = gz*(bes(1)-bes(3))/2. besr(2,1) = 0. besr(2,2) = gz*(bes(0)+bes(2))/2. besr(2,3) = gz*(bes(1)+bes(3))/4. besr(3,1) = gz*bes(0) besr(3,2) = gz*bes(1) besr(3,3) = gz*bes(2) besrd(1,1) = besr(2,2) besrd(1,2) = besr(2,3) besrd(2,1) = besr(1,2) besrd(2,2) = besr(1,3) end c/////////////////////////////////////////////////////////////////////// subroutine mYdtot(z,j,m,d,gz, nti,ti, Ydtot) c--------------------------------------------------------- c "d" means dash in the text c including subroutine: matse2 c matVd, matFd, msdeld (j=1, zd) c matVTd, msdeTd, matYsd (j=2) c--------------------------------------------------------- implicit double precision(a-h,o-z) dimension ti(nti) dimension Ydtot(2,2,nti), W1(2,2),W0(2,2), b(0:1) dimension Vd(2,2),Vds(2,2),VTd(2,2) dimension Fd(2,2),Dd(2,2),Dds(2,2),Ysd(2,2) dimension pW1(2,2),pW0(2,2) common /C1/ HL common /D2/ tau coef1 = exp(-abs(d-z)*gz) do 5 it=1,nti do 5 k=1,2 do 5 i=1,2 Ydtot(i,k,it) = 0. 5 continue c: (1)-(i) if (j.eq.1) then if (z.le.d) then call matVd(m,d,gz, Vd, 0) call matFd(z,gz, Fd) call matse2(Fd,Vd, W1) if (m.eq.1) then call matVd(m,d,gz, Vds, 1) call matse2(Fd,Vds, W0) else if (m.eq.2) then do 10 k=1,2 do 10 i=1,2 W0(i,k) = 0. 10 continue end if call msdeld(gz, b) c: (1)-(ii) else if (z.gt.d) then call matVTd(m,d,gz, VTd, 0) sw1=z-HL call matFd(sw1,gz, Fd) call matDd(Dd, -1, 0) call matDd(Dds, -1, 1) call matse2(Fd,Dd, pW1) call matse2(pW1,VTd, W1) call matse2(Fd,Dds, pW0) call matse2(pW0,VTd, W0) call msdeTd(gz, b) endif do 60 it=1,nti do 60 k=1,2 do 60 i=1,2 Ydtot(i,k,it) = coef1* (W0(i,k)/b(0) + & (W1(i,k)/b(1) - W0(i,k)/b(0))*exp(-b(0)/b(1)*ti(it)/tau)) 60 continue c: (2) else if (j.eq.2) then call msdeTd(gz, b) call matVTd(m,d,gz, W1, 0) if (m.eq.1) then do 110 it=1,nti do 110 k=1,2 do 110 i=1,2 Ydtot(i,k,it) = coef1* (W1(i,k)/b(0) + & (W1(i,k)/b(1) - W1(i,k)/b(0))*exp(-b(0)/b(1)*ti(it)/tau)) 110 continue else if (m.eq.2) then call matVTd(m,d,gz, W0, 1) coef2 = exp(-(z+d-2*HL)*gz) c index = sign(1,z-d) index = idnint(dsign(1.d0,z-d)) call matYsd(Ysd, index) do 120 it=1,nti do 120 k=1,2 do 120 i=1,2 Ydtot(i,k,it) = coef1*Ysd(i,k) + coef2* (W0(i,k)/b(0) + & (W1(i,k)/b(1) - W0(i,k)/b(0))*exp(-b(0)/b(1)*ti(it)/tau)) 120 continue endif endif end c----------------------------------------------------------------------- subroutine matVd(m,d,gz, Vd, is) c----------------------------------------------------------- c including subroutine: matDd, matFd, matYsd c----------------------------------------------------------- implicit double precision(a-h,o-z) dimension Vd(2,2), Qd(2,2), Dd(2,2),Fd(2,2) common /C1/ HL if (m.eq.2) then call matYsd(Qd, -1) else if (m.eq.1) then call matDd(Dd, 1, is) w1=HL-d call matFd(w1,gz, Fd) do 10 k=1,2 do 10 i=1,2 Qd(i,k) = 0. do 20 j1=1,2 Qd(i,k) = Qd(i,k) - 2*Dd(i,j1)*Fd(j1,k) 20 continue 10 continue endif do 50 k=1,2 Vd(1,k) = Qd(1,k)+Qd(2,k) Vd(2,k) = 0. 50 continue end c----------------------------------------------------------------------- subroutine matVTd(m,d,gz, VTd, is) c----------------------------------------------------------- c including subroutine: matFd, matPTd, matYsd, matse2 c----------------------------------------------------------- implicit double precision(a-h,o-z) dimension VTd(2,2), QTd(2,2), Fd(2,2),PTd(2,2),Ysd(2,2) common /C1/ HL if (m.eq.1) then call matFd(-d,gz, Fd) do 10 k=1,2 do 10 i=1,2 QTd(i,k) = -2*Fd(i,k) 10 continue else if (m.eq.2) then call matPTd(gz, PTd, is) call matYsd(Ysd, -1) call matse2(PTd,Ysd, QTd) else write(*,*) 'matQTd error' endif do 50 k=1,2 VTd(1,k) = -QTd(2,k) VTd(2,k) = -VTd(1,k) 50 continue end c----------------------------------------------------------------------- subroutine msdeld(gz, sdeld) c----------------------------------------------------------- c msdeld: make small delta dash c including subroutine: matPd, matQd c----------------------------------------------------------- implicit double precision(a-h,o-z) dimension sdeld(0:1), Pd(2,2),Pds(2,2) call matPd(gz, Pd, 0) call matPd(gz, Pds, 1) sdeld(1) = Pd(1,1)+Pd(2,1) sdeld(0) = Pds(1,1)+Pds(2,1) end c----------------------------------------------------------------------- subroutine msdeTd(gz, sdelTd) c----------------------------------------------------------- c msdeTd: make small delta tilde dash c including subroutine: matPTd, matQTd, matYsd c----------------------------------------------------------- implicit double precision(a-h,o-z) dimension sdelTd(0:1),PTd(2,2),PTds(2,2) call matPTd(gz, PTd, 0) call matPTd(gz, PTds, 1) sdelTd(1) = PTd(2,1)-PTd(2,2) sdelTd(0) = PTds(2,1)-PTds(2,2) end c----------------------------------------------------------------------- subroutine matPd(gz, Pd, is) c----------------------------------------------------------- c including subroutine: matDd, matFd c----------------------------------------------------------- implicit double precision(a-h,o-z) dimension Pd(2,2) dimension Fd1(2,2),Dd1(2,2) common /C1/ HL call matDd(Dd1, 1, is) call matFd(HL,gz, Fd1) call matse2(Dd1,Fd1, Pd) end c----------------------------------------------------------------------- subroutine matPTd(gz, PTd, is) c----------------------------------------------------------- c "T" corresponds to tilde (or Bar) in the text c including subroutine: matDd, matFd c----------------------------------------------------------- implicit double precision(a-h,o-z) dimension PTd(2,2) dimension FTd1(2,2),Dd1(2,2) common /C1/ HL call matDd(Dd1, -1, is) call matFd(-HL,gz, FTd1) call matse2(FTd1,Dd1, PTd) end c----------------------------------------------------------------------- subroutine matYsd(Ysd, index) implicit double precision(a-h,o-z) dimension Ysd(2,2) if (index.eq.-1) then Ysd(1,1) = -1. Ysd(1,2) = -1. Ysd(2,1) = -1. Ysd(2,2) = -1. else if (index.eq.1) then Ysd(1,1) = 1. Ysd(1,2) = -1. Ysd(2,1) = -1. Ysd(2,2) = 1. else write(*,*) 'matYsd error' stop endif end c----------------------------------------------------------------------- subroutine matFd(z,gz, Fd) implicit double precision(a-h,o-z) dimension Fd(2,2) zgz = gz*z Ch=(1+exp(-2*abs(zgz)))/2.0 Sh=(1-exp(-2*abs(zgz)))*idnint(dsign(1.d0,z))/2.0 Fd(1,1) = Ch Fd(1,2) = Sh Fd(2,1) = Sh Fd(2,2) = Ch end c----------------------------------------------------------------------- subroutine matDd(Dd, index, is) c----------------------------------------------------------- c index=-1 : inverse matrix c----------------------------------------------------------- implicit double precision(a-h,o-z) dimension Dd(2,2) common /B2/ um(2) do 10 i=1,2 do 10 ii=1,2 Dd(i,ii) = 0. 10 continue if (index.eq.1) then Dd(2,2) = um(1)/um(2) if(is.eq.0) Dd(1,1) = 1.0 else if (index.eq.-1) then Dd(1,1) = 1.0 if(is.eq.0) Dd(2,2) = um(2)/um(1) else write(*,*) 'D index error' stop end if if ((is.ne.0) .and. (is.ne.1)) then write(*,*) 'D is error' stop end if end c----------------------------------------------------------------------- subroutine matse2(A,B, C) c----------------------------------------------------------- c matse2 : matrix "seki" for 2*2 matrices c C = A*B c----------------------------------------------------------- implicit double precision(a-h,o-z) dimension C(2,2), A(2,2),B(2,2) do 10 i=1,2 do 10 j=1,2 C(i,j) = 0. do 10 j1=1,2 C(i,j) = C(i,j) + A(i,j1)*B(j1,j) 10 continue end c/////////////////////////////////////////////////////////////////////// subroutine mYtot(z,j,m,d,gz, nti,ti, Ytot) c--------------------------------------------------------- c including subroutine: Laplac, ppoly, divide c matVs, matF, matG (j=1, zd) c matVsT, matEzv, matYsv (j=2) c--------------------------------------------------------- implicit double precision(a-h,o-z) dimension ti(nti) dimension Ytot(4,3,nti), W(4,3,0:4), dVs(4,3,0:3),VsT(4,3,0:4) dimension F(4,4),G(4,4) dimension Dv(4,4,0:1),Env(4,4,0:1) dimension DE(4,4,0:2),DEV(4,3,0:6),dDEV(4,3,0:5) dimension Ysv(4,3,0:1),Ezv(4,4,0:1), EzVs(4,3,0:5),dEzVs(4,3,0:4) common /C1/ HL common /D1/ ak common /D2/ tau coef1 = exp(-abs(d-z)*gz) do 5 ns=0,4 do 5 k=1,3 do 5 i=1,4 W(i,k,ns) = 0. 5 continue c: (1)-(i) if (j.eq.1) then if (z.le.d) then call matVs(m,d,gz, dVs) call matF(z,0.0d0,j,gz, F) call matG(gz, G, 1) if (m.eq.2) then in = 1 else in = 0 end if do 10 ns=0,3 do 10 k=1,3 do 10 i=1,4 do 10 j1=1,4 do 10 j2=1,4 W(i,k,ns+in) = W(i,k,ns+in) & + F(i,j1)*G(j1,j2)*dVs(j2,k,ns) 10 continue call Laplac(W,m,gz, nti,ti, Ytot, 1) c: (1)-(ii) else if (z.gt.d) then call matVsT(m,d,gz, VsT) call matF(z,HL,j,gz, F) call matDv(Dv, -1) call matEnv(Env, 1) call ppoly(Dv,4,4,1,Env,4,4,1, DE) call ppoly(DE,4,4,2,VsT,4,3,4, DEV) call divide(DEV,6, dDEV) do 40 ns=0,3 do 40 k=1,3 do 40 i=1,4 do 40 j1=1,4 W(i,k,ns) = W(i,k,ns) + F(i,j1)*dDEV(j1,k,ns) 40 continue call Laplac(W,m,gz, nti,ti, Ytot, 2) endif do 60 it=1,nti do 60 k=1,3 do 60 i=1,4 Ytot(i,k,it) = Ytot(i,k,it)*coef1 60 continue c: (2)-(i) else if (j.eq.2) then call matVsT(m,d,gz, VsT) call matEzv(z,gz, Ezv) call ppoly(Ezv,4,4,1,VsT,4,3,4, EzVs) call divide(EzVs,5, dEzVs) if (m.eq.1) then do 100 k=1,3 do 100 i=1,4 W(i,k,3) = dEzVs(i,k,2) W(i,k,2) = dEzVs(i,k,2) + dEzVs(i,k,1) W(i,k,1) = dEzVs(i,k,1) + dEzVs(i,k,0) W(i,k,0) = dEzVs(i,k,0) 100 continue call Laplac(W,m,gz, nti,ti, Ytot, 2) do 110 it=1,nti do 110 k=1,3 do 110 i=1,4 Ytot(i,k,it) = Ytot(i,k,it)*coef1 110 continue c: (2)-(ii) else if (m.eq.2) then coef2 = exp(-(z+d-2*HL)*gz) call matYsv(z,d,gz, Ysv) call Laplac(dEzVs,m,gz, nti,ti, Ytot, 2) do 120 it=1,nti do 120 k=1,3 do 120 i=1,4 Ytot(i,k,it) = Ytot(i,k,it)*coef2 + coef1* (Ysv(i,k,0)/ak+ & (Ysv(i,k,1) - Ysv(i,k,0)/ak)*exp(-ak*ti(it)/tau)) 120 continue endif endif end c----------------------------------------------------------------------- subroutine Laplac(W,m,gz, nti,ti, Y, index) c---------------------------------------------------- c including subroutine: cardan, msdel, msdelT c index = 1 : down-going c 2 : up-going c---------------------------------------------------- implicit double precision(a-h,o-z) dimension ti(nti) dimension Y(4,3,nti), W(4,3,0:4),b(0:4),db(0:3),del(0:4) dimension Ye(4,3),zet(4) dimension den(4),up(4,3),C(4,3,7) common /D1/ ak common /D2/ tau if (index.eq.1) then call msdel(gz, b,db) else if (index.eq.2) then call msdelT(gz, b,db) else write(*,*) 'Laplac error' end if call cardan(db, zet) if (m.eq.1) then do 200 ns=0,3 del(ns) = db(ns) 200 continue zet(4) = 0. del(4) = 0. ly = 3 else if (m.eq.2) then do 210 ns=0,4 del(ns) = b(ns) 210 continue zet(4) = -ak ly = 4 end if do 10 i=1,4 do 10 k=1,3 Ye(i,k) = W(i,k,ly)/del(ly) c --- Ye: Elastic Solution --- 10 continue do 20 iv=1,ly den(iv) = 1. do 30 kv=1,ly if(kv.ne.iv) den(iv) = den(iv)*(zet(iv)-zet(kv)) 30 continue do 50 i=1,4 do 50 k=1,3 up(i,k) = 0. do 60 kw=0,ly-1 up(i,k) = up(i,k) + & (W(i,k,kw) - Ye(i,k)*del(kw))/del(ly)*((zet(iv))**kw) 60 continue C(i,k,iv) = up(i,k)/den(iv) c write(*,*) 'ly,gz, C= ',ly,real(gz),real(C(i,k,iv)) 50 continue 20 continue do 100 it=1,nti do 100 i=1,4 do 100 k=1,3 w1 = 0. do 110 iv=1,ly w1 = w1+ (1-exp(zet(iv)*ti(it)/tau))*C(i,k,iv)/zet(iv) 110 continue Y(i,k,it) = Ye(i,k) - w1 c write(11,*) 'Y,Ye ',i,k,real(Y(i,k)),real(Ye(i,k)) 100 continue end c----------------------------------------------------------------------- subroutine cardan(a, x) c---------------------------------------------------- c cardan: Cardano's fomula: solution of a 3 order equation: c a0 + a1*x + a2*xx + a3*xxx = 0 c input: a(i) c output (roots): x(i) c---------------------------------------------------- implicit double precision(a-h,o-z) parameter (pi=3.14159265358979, eps=1.0d-12) dimension x(3), a(0:3) w1= a(2)/3/a(3) p = a(1)/3/a(3)-w1*w1 q = (a(0)-a(1)*w1)/a(3) + 2*w1**3 s = -(q*q + 4*p**3) if (abs(s).lt.eps) then if (q.ge.0) then x(1) = -2*(q/2.)**(1.0/3.0) -w1 x(2) = (q/2.)**(1.0/3.0) -w1 else x(1) = 2*(-q/2.)**(1.0/3.0) -w1 x(2) = -(-q/2.)**(1.0/3.0) -w1 end if x(3) = x(2) write(*,*) 'Cardano; the same roots' else if (s.gt.0.) then r = dsqrt(q*q+s)/2 the= datan2(sqrt(s),-q) r3 = 2*r**(1.0/3.0) x(1) = r3*dcos(the/3.) -w1 x(2) = r3*dcos((the+2*pi)/3.) -w1 x(3) = r3*dcos((the+4*pi)/3.) -w1 else write(*,*) 'Cardan error; Complex solution' end if end c----------------------------------------------------------------------- subroutine divide(A,n, dA) c--------------------------------------------------------- c divide : divide a matrix A (4*3) by ak c dA = A/ak c--------------------------------------------------------- implicit double precision(a-h,o-z) dimension A(4,3,0:n),dA(4,3,0:n-1) common /D1/ ak do 100 k=1,3 do 100 i=1,4 dA(i,k,n-1) = A(i,k,n) dA(i,k,0) = A(i,k,0)/ak do 110 ns=n-2,1,-1 dA(i,k,ns) = A(i,k,ns+1) - ak*dA(i,k,ns+1) 110 continue 100 continue end c----------------------------------------------------------------------- c######## CHECK ### c c subroutine matY0(m,d,gz, Y0) c implicit double precision(a-h,o-z) c dimension Y0(4,3), Vs(4,3,0:4),sdel(0:4),dsdel(0:3) c c call matVs(m,d,gz, Vs) c call msdel(gz, sdel,dsdel) c do 10 k=1,3 c do 10 i=1,4 c Y0(i,k) = Vs(i,k,4)/sdel(4) c 10 continue c end cc--- c subroutine matAn(m,d,gz, An) c implicit double precision(a-h,o-z) c dimension An(4,3), VsT(4,3,0:4),sdelT(0:4),dsdelT(0:3) c c call matVsT(m,d,gz, VsT) c call msdelT(gz, sdelT,dsdelT) c do 10 k=1,3 c do 10 i=1,4 c if(m.eq.2) then c An(i,k) = VsT(i,k,4)/sdelT(4) c else if (m.eq.1) then c An(i,k) = VsT(i,k,2)/sdelT(4) c endif c 10 continue c end c c####### CHECK ### c----------------------------------------------------------------------- subroutine matVs(m,d,gz, dVs) c----------------------------------------------------------- c including subroutine: scaP, vecQ, ppolyv c----------------------------------------------------------- implicit double precision(a-h,o-z) dimension Vs(4,3,0:4),dVs(4,3,0:3) dimension Pa(0:2),Pb(0:2),Pc(0:2),Pd(0:2) dimension Qa(3,0:2),Qb(3,0:2) dimension PQ1(3,0:4),PQ2(3,0:4),PQ3(3,0:4),PQ4(3,0:4) call scaP(gz, Pa,Pb,Pc,Pd) call vecQ(m,d,gz, Qa,Qb) call ppolyv(Pd,2,Qa,2, PQ1) call ppolyv(Pb,2,Qb,2, PQ2) call ppolyv(Pc,2,Qa,2, PQ3) call ppolyv(Pa,2,Qb,2, PQ4) do 50 ns=0,4 do 50 k=1,3 Vs(1,k,ns) = PQ1(k,ns) - PQ2(k,ns) Vs(2,k,ns) = -PQ3(k,ns) + PQ4(k,ns) Vs(3,k,ns) = 0. Vs(4,k,ns) = 0. 50 continue call divide(Vs,4, dVs) end c----------------------------------------------------------------------- subroutine matVsT(m,d,gz, VsT) c--------------------------------------------------------- c including subroutine: scaPT, vecQT, ppolyv c--------------------------------------------------------- implicit double precision(a-h,o-z) dimension VsT(4,3,0:4) dimension PTa(0:2),PTb(0:2),PTc(0:2),PTd(0:2) dimension QTa(3,0:2),QTb(3,0:2) dimension PQ1(3,0:4),PQ2(3,0:4),PQ3(3,0:4),PQ4(3,0:4) call scaPT(gz, PTa,PTb,PTc,PTd) call vecQT(m,d,gz, QTa,QTb) call ppolyv(PTd,2,QTa,2, PQ1) call ppolyv(PTb,2,QTb,2, PQ2) call ppolyv(PTc,2,QTa,2, PQ3) call ppolyv(PTa,2,QTb,2, PQ4) do 10 ns=0,4 do 10 k=1,3 VsT(1,k,ns) = -PQ1(k,ns) + PQ2(k,ns) VsT(2,k,ns) = PQ3(k,ns) - PQ4(k,ns) VsT(3,k,ns) = -VsT(1,k,ns) VsT(4,k,ns) = -VsT(2,k,ns) 10 continue end c----------------------------------------------------------------------- subroutine msdel(gz, sdel,dsdel) c---------------------------------------------------- c msdel : make small-delta and small reverse-prime delta c including subroutine: scaP, ppoly0 c---------------------------------------------------- implicit double precision(a-h,o-z) dimension sdel(0:4),dsdel(0:3) dimension Pa(0:2),Pb(0:2),Pc(0:2),Pd(0:2), PE(0:4),PF(0:4) common /D1/ ak call scaP(gz, Pa,Pb,Pc,Pd) call ppoly0(Pa,2,Pd,2, PE) call ppoly0(Pb,2,Pc,2, PF) do 10 ns=0,4 sdel(ns) = PE(ns) - PF(ns) 10 continue c devide "sdel" by "S+ak" dsdel(3) = sdel(4) dsdel(2) = sdel(3) - ak*dsdel(3) dsdel(1) = sdel(2) - ak*dsdel(2) dsdel(0) = sdel(0)/ak end c----------------------------------------------------------------------- subroutine msdelT(gz, sdelT,dsdelT) c---------------------------------------------------- c msdelT : make small-delta-Tilde and reverse-prime delta Tilde c including subroutine: scaPT, ppoly0 c---------------------------------------------------- implicit double precision(a-h,o-z) dimension sdelT(0:4),dsdelT(0:3) dimension PTa(0:2),PTb(0:2),PTc(0:2),PTd(0:2), PTE(0:4),PTF(0:4) common /D1/ ak call scaPT(gz, PTa,PTb,PTc,PTd) call ppoly0(PTa,2,PTd,2, PTE) call ppoly0(PTb,2,PTc,2, PTF) do 10 ns=0,4 sdelT(ns) = PTE(ns) - PTF(ns) 10 continue dsdelT(3) = sdelT(4) dsdelT(2) = sdelT(3) - ak*dsdelT(3) dsdelT(1) = sdelT(2) - ak*dsdelT(2) dsdelT(0) = sdelT(0)/ak end c----------------------------------------------------------------------- subroutine scaP(gz, Pa,Pb,Pc,Pd) c---------------------------------------------------- c including subroutine: matEnv, matDv, matF, matG, ppoly c---------------------------------------------------- implicit double precision(a-h,o-z) dimension P(4,4,0:2), Pa(0:2),Pb(0:2),Pc(0:2),Pd(0:2) dimension Env(4,4,0:1),Dv(4,4,0:1), F1(4,4),G(4,4) dimension ED(4,4,0:2) common /C1/ HL call matEnv(Env, -1) call matDv(Dv, 1) call ppoly(Env,4,4,1,Dv,4,4,1, ED) call matF(HL,0.0d0,1,gz, F1) call matG(gz, G, 1) do 10 ns=0,2 do 10 i=1,4 do 10 ii=1,4 P(i,ii,ns) = 0.0 do 10 j4=1,4 do 10 j5=1,4 P(i,ii,ns) = P(i,ii,ns)+ ED(i,j4,ns)*F1(j4,j5)*G(j5,ii) 10 continue do 50 ns=0,2 Pa(ns) = P(1,1,ns) + P(3,1,ns) Pb(ns) = P(1,2,ns) + P(3,2,ns) Pc(ns) = P(2,1,ns) + P(4,1,ns) Pd(ns) = P(2,2,ns) + P(4,2,ns) 50 continue end c----------------------------------------------------------------------- subroutine vecQ(m,d,gz, Qa,Qb) c---------------------------------------------------- c including subroutine: matEnv, matDv, matF, matDel, matYsv, ppoly c---------------------------------------------------- implicit double precision(a-h,o-z) dimension Q(4,3,0:2), Qa(3,0:2),Qb(3,0:2) dimension Env(4,4,0:1),Dv(4,4,0:1), ED(4,4,0:2) dimension F1(4,4), Del(4,3),Ysv(4,3,0:1) common /C1/ HL call matEnv(Env, -1) if(m.eq.1) then call matDv(Dv, 1) call ppoly(Env,4,4,1,Dv,4,4,1, ED) call matF(HL,d,1,gz, F1) call matDel(Del) do 10 ns=0,2 do 10 k=1,3 do 10 i=1,4 Q(i,k,ns) = 0.0 do 10 j2=1,4 do 10 j3=1,4 Q(i,k,ns) = Q(i,k,ns) + ED(i,j2,ns)*F1(j2,j3)*Del(j3,k) 10 continue else if(m.eq.2) then call matYsv(HL,d,gz, Ysv) call ppoly(Env,4,4,1,Ysv,4,3,1, Q) else write(*,*) 'matQ error' stop endif do 30 ns=0,2 do 30 k=1,3 Qa(k,ns) = Q(1,k,ns) + Q(3,k,ns) Qb(k,ns) = Q(2,k,ns) + Q(4,k,ns) 30 continue end c----------------------------------------------------------------------- subroutine scaPT(gz, PTa,PTb,PTc,PTd) c---------------------------------------------------- c "T" means tilde in the text c including subroutine: matEnv, matDv, matF, matG, ppoly c---------------------------------------------------- implicit double precision(a-h,o-z) dimension PT(4,4,0:2), PTa(0:2),PTb(0:2),PTc(0:2),PTd(0:2) dimension G(4,4),FT1(4,4) dimension Dv(4,4,0:1),Env(4,4,0:1), DE(4,4,0:2) common /C1/ HL call matG(gz, G, -1) call matF(0.0d0,HL,1,gz, FT1) call matDv(Dv, -1) call matEnv(Env, 1) call ppoly(Dv,4,4,1,Env,4,4,1, DE) do 10 ns=0,2 do 10 i=1,4 do 10 ii=1,4 PT(i,ii,ns) = 0.0 do 10 j1=1,4 do 10 j2=1,4 PT(i,ii,ns) = PT(i,ii,ns)+ G(i,j1)*FT1(j1,j2)*DE(j2,ii,ns) 10 continue do 50 ns=0,2 PTa(ns) = PT(3,1,ns) - PT(3,3,ns) PTb(ns) = PT(3,2,ns) - PT(3,4,ns) PTc(ns) = PT(4,1,ns) - PT(4,3,ns) PTd(ns) = PT(4,2,ns) - PT(4,4,ns) 50 continue end c----------------------------------------------------------------------- subroutine vecQT(m,d,gz, QTa,QTb) c---------------------------------------------------- c including subroutine: matG, matF, matDv, matYsv, matDel, ppoly c---------------------------------------------------- implicit double precision(a-h,o-z) dimension QT(4,3,0:2), QTa(3,0:2),QTb(3,0:2) dimension G(4,4),FT1(4,4),Del(4,3) dimension Dv(4,4,0:1),Ysv(4,3,0:1), DY(4,3,0:2) common /C1/ HL do 1 ns=0,2 do 1 k=1,3 do 1 i=1,4 QT(i,k,ns) = 0. 1 continue call matG(gz, G, -1) if (m.eq.1) then call matF(0.0d0,d,1,gz, FT1) call matDel(Del) do 5 k=1,3 do 5 i=1,4 QT(i,k,0) = 0. do 5 j1=1,4 do 5 j2=1,4 QT(i,k,0) = QT(i,k,0) + G(i,j1)*FT1(j1,j2)*Del(j2,k) 5 continue else if (m.eq.2) then call matF(0.0d0,HL,1,gz, FT1) call matDv(Dv, -1) call matYsv(HL,d,gz, Ysv) call ppoly(Dv,4,4,1,Ysv,4,3,1, DY) do 10 ns=0,2 do 10 k=1,3 do 10 i=1,4 QT(i,k,ns) = 0. do 15 j1=1,4 do 15 j2=1,4 QT(i,k,ns) = QT(i,k,ns) + G(i,j1)*FT1(j1,j2)*DY(j2,k,ns) 15 continue 10 continue endif do 50 ns=0,2 do 50 k=1,3 QTa(k,ns) = QT(3,k,ns) QTb(k,ns) = QT(4,k,ns) 50 continue end c----------------------------------------------------------------------- subroutine matYsv(z,d,gz, Ysv) c---------------------------------------------------- c including subroutine: matYs, makev c---------------------------------------------------- implicit double precision(a-h,o-z) dimension Ysv(4,3,0:1), Ys(4,3), Yss(4,3) call matYs(z,d,gz, Ys, 0) call matYs(z,d,gz, Yss,1) call makev(Ys,Yss,4,3, Ysv) end c----------------------------------------------------------------------- subroutine matYs(z,d,gz, Ys, is) implicit double precision(a-h,o-z) dimension Ys(4,3) common /B1/ gam(2) ga = gam(2) if (is.eq.1) ga = 1.0 dsg = gz*(z-d) if (z.lt.d) then Ys(1,1) = 2 + ga*(1+3*dsg) Ys(2,1) = -4 + ga*(4+3*dsg) Ys(3,1) = -1 + ga*(4+3*dsg) Ys(4,1) = -1 + ga*(1+3*dsg) Ys(1,2) = 1 + ga*dsg Ys(2,2) = -1 + ga*(1+dsg) Ys(3,2) = ga*(1+dsg) Ys(4,2) = ga*dsg Ys(1,3) = 2 + ga*(-1+dsg) Ys(2,3) = ga*dsg Ys(3,3) = 1 + ga*dsg Ys(4,3) = 1 + ga*(-1+dsg) else Ys(1,1) = 2 + ga*(1-3*dsg) Ys(2,1) = 4 + ga*(-4+3*dsg) Ys(3,1) = 1 + ga*(-4+3*dsg) Ys(4,1) = -1 + ga*(1-3*dsg) Ys(1,2) = -1 + ga*dsg Ys(2,2) = -1 + ga*(1-dsg) Ys(3,2) = ga*(1-dsg) Ys(4,2) = ga*dsg Ys(1,3) = 2 + ga*(-1-dsg) Ys(2,3) = ga*dsg Ys(3,3) = -1 + ga*dsg Ys(4,3) = 1 + ga*(-1-dsg) endif end c----------------------------------------------------------------------- subroutine matDel(Del) implicit double precision(a-h,o-z) dimension Del(4,3) common /B1/ gam(2) ga = gam(1) do 10 k=1,3 do 10 i=1,4 Del(i,k) = 0. 10 continue Del(2,1) = -8 + ga*8 Del(3,1) = -2 + ga*8 Del(1,2) = 2. Del(3,3) = 2. end c----------------------------------------------------------------------- subroutine matF(z1,z2,j,gz, F) c---------------------------------------------------- c matF : Fj(z)=F(j;z1-z2) c---------------------------------------------------- implicit double precision(a-h,o-z) dimension F(4,4) common /B1/ gam(2) ga = gam(j) z = z1-z2 zgz = gz*z Ch=(1+exp(-2*abs(zgz)))/2.0 Sh=(1-exp(-2*abs(zgz)))*idnint(dsign(1.d0,z))/2.0 F(1,1) = Ch + ga* zgz*Sh F(1,2) = -Sh + ga*( Sh-zgz*Ch) F(1,3) = 2*Sh + ga*(-Sh+zgz*Ch) F(1,4) = - ga* zgz*Sh F(2,1) = -Sh + ga*(zgz*Ch+Sh) F(2,2) = Ch - ga* zgz*Sh F(2,3) = ga* zgz*Sh F(2,4) = 2*Sh - ga*(zgz*Ch+Sh) F(3,1) = ga*(zgz*Ch+Sh) F(3,2) = - ga* zgz*Sh F(3,3) = Ch + ga* zgz*Sh F(3,4) = Sh - ga*(zgz*Ch+Sh) F(4,1) = ga* zgz*Sh F(4,2) = + ga*( Sh-zgz*Ch) F(4,3) = Sh + ga*(-Sh+zgz*Ch) F(4,4) = Ch - ga* zgz*Sh end c----------------------------------------------------------------------- subroutine matEzv(z,gz, Ezv) c---------------------------------------------------- c including subroutine: matEnz, makev c---------------------------------------------------- implicit double precision(a-h,o-z) dimension Ezv(4,4,0:1), Ez(4,4),Ezs(4,4) call matEnz(z,gz, Ez, 0) call matEnz(z,gz, Ezs,1) call makev(Ez,Ezs,4,4, Ezv) end c----------------------------------------------------------------------- subroutine matEnz(z,gz, Enz, is) c---------------------------------------------------- c matEnz : E(n+1;z)=Enz(z-Hn); Hn==HL(nt-1) c---------------------------------------------------- implicit double precision(a-h,o-z) dimension Enz(4,4) common /B1/ gam(2) common /C1/ HL ga = gam(2) if (is.eq.1) ga = 1.0 zn = z-HL ggz = ga*zn*gz Enz(1,1) = 1. Enz(1,2) = 0. Enz(1,3) = 0. Enz(1,4) = ggz Enz(2,1) = 0. Enz(2,2) = 2-ga + ggz Enz(2,3) = 1. Enz(2,4) = 0. Enz(3,1) = 0. Enz(3,2) = 1-ga + ggz Enz(3,3) = 1. Enz(3,4) = 0. Enz(4,1) = 1. Enz(4,2) = 0. Enz(4,3) = 0. Enz(4,4) = 1 + ggz end c----------------------------------------------------------------------- subroutine matEnv(Env, index) c---------------------------------------------------- c including subroutine: matEn, makev c---------------------------------------------------- implicit double precision(a-h,o-z) dimension Env(4,4,0:1), En(4,4),Es(4,4) call matEn(En, index, 0) call matEn(Es, index, 1) call makev(En,Es,4,4, Env) end c----------------------------------------------------------------------- subroutine makev(A,As,nr,nc, Av) implicit double precision(a-h,o-z) dimension Av(nr,nc,0:1), A(nr,nc),As(nr,nc) common /D1/ ak do 10 i=1,nr do 10 ii=1,nc Av(i,ii,1) = A(i,ii) Av(i,ii,0) = ak*As(i,ii) 10 continue end c----------------------------------------------------------------------- subroutine matEn(En, index, is) c---------------------------------------------------- c E(n+1) = Enz(n+1;z=0) (in this program n+1=nt, E(n+1)=En) c index: -1 means inverse matrix c including subroutine: matEnz c---------------------------------------------------- implicit double precision(a-h,o-z) dimension En(4,4) common /B1/ gam(2) common /C1/ HL ga = gam(2) if (is.eq.1) ga = 1.0 call matEnz(HL,0.0d0, En, is) if (index.eq.-1) then En(2,2) = 1. En(2,3) = -1. En(3,2) = -1+ga En(3,3) = 2-ga En(4,1) = -1. end if end c----------------------------------------------------------------------- subroutine matDv(Dv, index) c---------------------------------------------------- c including subroutine: matDL c---------------------------------------------------- implicit double precision(a-h,o-z) dimension Dv(4,4,0:1), DL(4,4),Ds(4,4) call matDL(DL, index, 0) call matDL(Ds, index, 1) do 10 i=1,4 do 10 ii=1,4 Dv(i,ii,1) = DL(i,ii) Dv(i,ii,0) = Ds(i,ii) 10 continue end c----------------------------------------------------------------------- subroutine matDL(DL, index, is) c---------------------------------------------------- c index: -1 means inverse matrix c---------------------------------------------------- implicit double precision(a-h,o-z) dimension DL(4,4) common /B2/ um(2) do 10 i=1,4 do 10 ii=1,4 DL(i,ii) = 0.0 10 continue if (index.eq.1) then DL(3,3) = um(1)/um(2) if (is.eq.0) DL(1,1) = 1.0 else if (index.eq.-1) then DL(1,1) = 1.0 if (is.eq.0) DL(3,3) = um(2)/um(1) else write(*,*) 'D index error' stop end if if ((is.ne.0) .and. (is.ne.1)) then write(*,*) 'D is error' stop end if DL(2,2) = DL(1,1) DL(4,4) = DL(3,3) end c----------------------------------------------------------------------- subroutine matG(gz, G, index) implicit double precision(a-h,o-z) dimension G(4,4) common /B2/ um(2) common /B3/ rou(2) common /B4/ gra G(1,1) = 1. G(1,2) = 0. G(1,3) = 0. G(1,4) = 0. G(2,1) = 0. G(2,2) = 1. G(2,3) = 0. G(2,4) = 0. G(3,1) = 0. G(3,2) = 0. G(3,3) = 1. G(3,4) = 0. G(4,1) = 0. G(4,2) = rou(1)*gra/2./um(1)/gz G(4,3) = 0. G(4,4) = 1. if (index.eq.-1) G(4,2)=-G(4,2) end c----------------------------------------------------------------------- subroutine ppoly0(A,n1,B,n2, C) c---------------------------------------------------- c pro_poly : Products of Polynominals for numbers c A : number; A(0) + A(1)*x + ... + A(n1)*x**n1 c B : number; B(0) + B(1)*x + ... + B(n1)*x**n1 c---------------------------------------------------- implicit double precision(a-h,o-z) dimension C(0:n1+n2),A(0:n1),B(0:n2) do 50 ip=0,n1+n2 C(ip) = 0. call rangep(ip,n1,n2, inti,iend) do 100 ik=inti,iend C(ip) = C(ip) + A(ip-ik)*B(ik) 100 continue 50 continue end c----------------------------------------------------------------------- subroutine ppolyv(A,n1,B,n2, C) c---------------------------------------------------- c pro_poly : Products of Polynominals for (number)*(vector) c A : number c B : vector c---------------------------------------------------- implicit double precision(a-h,o-z) dimension C(3,0:n1+n2),A(0:n1),B(3,0:n2) do 50 ip=0,n1+n2 call rangep(ip,n1,n2, inti,iend) do 100 k=1,3 C(k,ip) = 0. do 100 ik=inti,iend C(k,ip) = C(k,ip) + A(ip-ik)*B(k,ik) 100 continue 50 continue end c----------------------------------------------------------------------- subroutine ppoly(A,nAr,nAc,nsA,B,nBr,nBc,nsB, C) c---------------------------------------------------- c ppoly : Products of Polynominals for (matrix)*(matrix) c input: A, B c A : nAr*nAc matrix; A(0) + A(1)*x + ... + A(nsA)*x**n1 c B : nBr*nBc matrix; B(0) + B(1)*x + ... + B(nsB)*x**n1 c r: row, c: column c Caution!! A and B is not changeable c c output: C : (nsA+nsB) order Polynominal & (nAr*nBc) matrix c---------------------------------------------------- implicit double precision(a-h,o-z) dimension A(nAr,nAc,0:nsA),B(nBr,nBc,0:nsB) dimension C(nAr,nBc,0:nsA+nsB) if (nAc.ne.nBr) then write(*,*) "ppoly Error" stop endif c write(*,*) 'B=',((real(B(i,k,6)),i=1,4),k=1,3) do 10 ip=0,nsA+nsB do 10 i=1,nAr do 10 k=1,nBc C(i,k,ip) = 0. 10 continue c write(*,*) '2B=',((real(B(i,k,6)),i=1,4),k=1,3) do 50 ip=0,nsA+nsB call rangep(ip,nsA,nsB, inti,iend) do 100 ik=inti,iend do 100 i=1,nAr do 100 k=1,nBc do 100 j1=1,nAc C(i,k,ip) = C(i,k,ip) + A(i,j1,ip-ik)*B(j1,k,ik) 100 continue 50 continue end c----------------------------------------------------------------------- subroutine rangep(ip,n1,n2, inti,iend) implicit double precision(a-h,o-z) if (ip.le.n1) then inti = 0 else inti = ip-n1 endif if (ip.le.n2) then iend = ip else iend = n2 endif end c/////////////////////////////////////////////////////////////////////// FUNCTION bessj0(x) implicit double precision(a-h,o-z) DOUBLE PRECISION bessj0,x DOUBLE PRECISION ax,xx,z DOUBLE PRECISION p1,p2,p3,p4,p5,q1,q2,q3,q4,q5,r1,r2,r3,r4,r5,r6, *s1,s2,s3,s4,s5,s6,y SAVE p1,p2,p3,p4,p5,q1,q2,q3,q4,q5,r1,r2,r3,r4,r5,r6,s1,s2,s3,s4, *s5,s6 DATA p1,p2,p3,p4,p5/1.d0,-.1098628627d-2,.2734510407d-4, *-.2073370639d-5,.2093887211d-6/, q1,q2,q3,q4,q5/-.1562499995d-1, *.1430488765d-3,-.6911147651d-5,.7621095161d-6,-.934945152d-7/ DATA r1,r2,r3,r4,r5,r6/57568490574.d0,-13362590354.d0, *651619640.7d0,-11214424.18d0,77392.33017d0,-184.9052456d0/,s1,s2, *s3,s4,s5,s6/57568490411.d0,1029532985.d0,9494680.718d0, *59272.64853d0,267.8532712d0,1.d0/ if(abs(x).lt.8.)then y=x**2 bessj0=(r1+y*(r2+y*(r3+y*(r4+y*(r5+y*r6)))))/(s1+y*(s2+y*(s3+y* *(s4+y*(s5+y*s6))))) else ax=abs(x) z=8./ax y=z**2 xx=ax-.785398164 bessj0=sqrt(.636619772/ax)*(cos(xx)*(p1+y*(p2+y*(p3+y*(p4+y* *p5))))-z*sin(xx)*(q1+y*(q2+y*(q3+y*(q4+y*q5))))) endif return END C (C) Copr. 1986-92 Numerical Recipes Software '035a%10(3z1+91;*0. c----------------------------------------------------------------------- FUNCTION bessj1(x) implicit double precision(a-h,o-z) DOUBLE PRECISION bessj1,x DOUBLE PRECISION ax,xx,z DOUBLE PRECISION p1,p2,p3,p4,p5,q1,q2,q3,q4,q5,r1,r2,r3,r4,r5,r6, *s1,s2,s3,s4,s5,s6,y SAVE p1,p2,p3,p4,p5,q1,q2,q3,q4,q5,r1,r2,r3,r4,r5,r6,s1,s2,s3,s4, *s5,s6 DATA r1,r2,r3,r4,r5,r6/72362614232.d0,-7895059235.d0, *242396853.1d0,-2972611.439d0,15704.48260d0,-30.16036606d0/,s1,s2, *s3,s4,s5,s6/144725228442.d0,2300535178.d0,18583304.74d0, *99447.43394d0,376.9991397d0,1.d0/ DATA p1,p2,p3,p4,p5/1.d0,.183105d-2,-.3516396496d-4, *.2457520174d-5,-.240337019d-6/, q1,q2,q3,q4,q5/.04687499995d0, *-.2002690873d-3,.8449199096d-5,-.88228987d-6,.105787412d-6/ if(abs(x).lt.8.)then y=x**2 bessj1=x*(r1+y*(r2+y*(r3+y*(r4+y*(r5+y*r6)))))/(s1+y*(s2+y*(s3+ *y*(s4+y*(s5+y*s6))))) else ax=abs(x) z=8./ax y=z**2 xx=ax-2.356194491 bessj1=sqrt(.636619772/ax)*(cos(xx)*(p1+y*(p2+y*(p3+y*(p4+y* *p5))))-z*sin(xx)*(q1+y*(q2+y*(q3+y*(q4+y*q5)))))*dsign(1.d0,x) endif return END C (C) Copr. 1986-92 Numerical Recipes Software '035a%10(3z1+91;*0. c----------------------------------------------------------------------- FUNCTION bessj(n,x) implicit double precision(a-h,o-z) INTEGER n,IACC DOUBLE PRECISION bessj,x,BIGNO,BIGNI PARAMETER (IACC=40,BIGNO=1.d10,BIGNI=1.d-10) CU USES bessj0,bessj1 INTEGER j,jsum,m DOUBLE PRECISION ax,bj,bjm,bjp,sum,tox,bessj0,bessj1 if(n.lt.2)pause 'bad argument n in bessj' ax=abs(x) if(ax.eq.0.)then bessj=0. else if(ax.gt.float(n))then tox=2./ax bjm=bessj0(ax) bj=bessj1(ax) do 11 j=1,n-1 bjp=j*tox*bj-bjm bjm=bj bj=bjp 11 continue bessj=bj else tox=2./ax m=2*((n+int(sqrt(float(IACC*n))))/2) bessj=0. jsum=0 sum=0. bjp=0. bj=1. do 12 j=m,1,-1 bjm=j*tox*bj-bjp bjp=bj bj=bjm if(abs(bj).gt.BIGNO)then bj=bj*BIGNI bjp=bjp*BIGNI bessj=bessj*BIGNI sum=sum*BIGNI endif if(jsum.ne.0)sum=sum+bj jsum=1-jsum if(j.eq.n)bessj=bjp 12 continue sum=2.*sum-bj bessj=bessj/sum endif if(x.lt.0..and.mod(n,2).eq.1)bessj=-bessj return END C (C) Copr. 1986-92 Numerical Recipes Software '035a%10(3z1+91;*0.