duoceng 2.f90

fortran语言编写的

program duoceng

integer,parameter::d=24,cyd=126,cyy=188
real vprm(cyd),vsrm(cyd),denrm(cyd),bsbrm(cyd),vprmh(cyd),vprmc(cyd)
real bsbrmh(cyd),bsbrmc(cyd),denrmh(cyd),denrmc(cyd)
real,dimension(cyy)::vp,vs,tt0,h,vsc,vph,vpc,vsh,denh,denc,t

real cf(d,cyy),rad(d,cyy)
real,dimension(cyd)::wav,buff
real,dimension(200,500):: seis,ra
complex res
real:: l,ratio,samp,samptt,sampt,scaler,a1,pi2,w
real::f0,t0,wl,w0,lsyn,sss,c,hs,n,r0,r2,r4,r,vp1,vp2
real:: ccc,s,f,a0,b,m
integer::j,i,lw,lb,a
real ,parameter ::pi=3.1415926,x1=12.5,ox=12.5,crp=2.0
data vp(1:40 )/40*3000/
data vp(41:75)/35*3500/
data vp(76:126)/51*3000/




           bsb1=0.3
	       bsb2=0.3

	      bsb3=0.3
          samp  = 2.0
          samptt=samp/1000
          t(0)=0.
          h(0)=0.
   
   
    do j=1,cyd
	   if(j.ge.1.and.j.le.40)  then
            vs(j)=((0.5-bsb1)/(1-bsb1))*vp(j)
		 else if(j.ge.41.and.j.le.75)     then          !vs(i)
		    vs(j)=((0.5-bsb2)/(1-bsb2))*vp(j)
		 else
             vs(j)=((0.5-bsb3)/(1-bsb3))*vp(j)
		endif
	   h(j)=vp(j)*samptt+h(j-1)           !heigth  all
	   t(j)=samptt+t(j-1)                 !t(i)=h(j)/vp(j) 
	  
    enddo



   do i=1,d
         
     do j=1,cyd
         vpc(0)=0.0
         vsc(0)=0.0

         vpc(j)=samptt*vp(j)**2+vpc(j-1)
	     vprm(j)=sqrt(vpc(j)/t(j))
         vprmh(j)=(vprm(j+1)+vprm(j))/2.0
         vprmc(j)=vprm(j+1)-vprm(j)
      

         vsc(j)=samptt*vs(j)**2+vsc(j-1)
         vsrm(j)=sqrt(vsc(j)/t(j))

         bsbrm(j)=1-vprm(j)/(2.0*(vprm(j)**2-vsrm(j)**2))
         bsbrmh(j)=(bsbrm(j+1)+bsbrm(j))/2.0
         bsbrmc(j)=bsbrm(j+1)-bsbrm(j)
       

         denrm(j)=crp*(vprm(j)**2)**(1.0/4.0)
		 denrmh(j)=(denrm(j+1)+denrm(j))/2.0
         denrmc(j)=denrm(j+1)-denrm(j)

      if(j.eq.cyd) then
	       vprmh(j)=vprmh(j-1)
	       vprmc(j)=vprmc(j-1)
           bsbrmh(j)=bsbrmh(j-1) 
		   bsbrmc(j)=bsbrmc(j-1) 
           denrmh(j)=denrmh(j-1)
           denrmc(j)=denrmc(j-1)
	  endif
    !write(10,*) j*2,vprm(j),vsrm(j),denrm(j)
    enddo
enddo




                f0 = 35.0
	            t0 = 0.0
	            samp  = 2.0
	            wl = 126
				w0 = 0.0
	            samptt= samp/1000.0    
		        lw=wl/samp +0.5
			    c=1.0
			    sss=5.5                 	         
			 if( mod(lw*1.,2.).eq.0 ) then
                  lw = lw+1 
		   		  t0 = t0+samp*2.0   
             endif

			 print*,'lw====',lw

             call wavelet(f0,t0,w0,c,sss,samptt,wav,lw)	
	                lsyn = lb1+lw-1


         scaler=100.0
      do i=1,d
	            l=x1+ox*(i-1)
		  do j=1,cyd    
	            ra(i,j)=atan(l/(2.0*h(j)))
			    rad(i,j)=ra(i,j)*180/pi
			    hs=sin(rad(i,j))/vprm(j)
	         if (j.lt.40) then	         
                b=0.0
			 else			  
                b=(vprmc(j)/vprmh(j))/(vprmc(j)/vprmh(j)+denrmc(j)/denrmh(j))
		     endif
                a0=b-2.0*(1.0+b)*(1.0-2.0*bsbrmh(j))/(1.0-bsbrmh(j))
                m=a0*r0+bsbrmc(j)/((1.0-bsbrmh(j))**2)
                n=vprmc(j)/(vprmh(j)*2.0)
                r0=(vprmc(j)/vprmh(j)+denrmc(j)/denrmh(j))/2.0
                   
             call Rtcoef(vprmh(j),vprm(j),vprm(j+1),vprmc(j),hs,r0,m,n,res)
			     if(vprm(j+1).gt.vprm(j))  then
		            cf(i,j)=sqrt(real(res)**2+imag(res)**2)
				   else
				     cf(i,j)=-sqrt(real(res)**2+imag(res)**2)
				   endif
                   if(j.eq.cyd)          cf(i,j)=cf(i,j-1)
				   if(j.eq.1)            cf(i,j)=cf(i,j+1)
				   if(j.lt.1.or.j.gt.cyd)   cf(i,j)=0
               if(i.eq.1.) print*,'i===',i,j, ra(i,j),cf(i,j)
                
          enddo
              
             ! call fold(lb1,cf,lw,wav,lsyn,buff,scaler)
		         do j=1,lsyn
		             if(i.eq.5) print*,'i===',i,j,buff(j),cf(i,j)
		                 seis(i,j)=buff(j)
			!			 print*,seis(i,j)
	             enddo 
         
	  enddo	     
             

          open(unit=300,file='ra.dat',status='old')

          open(unit=100,file='reflection.dat',status='old')
!......................................................
        
           do j=1,cyy
                write(100,301)  j*2.0,(cf(i,j)+(i-1)*10,i=1,18)
                write(300,301)  j*2.0,(ra(i,j),i=1,18)
		   enddo
301            format(1x,f6.1,18f10.3)

        open(unit=300,file='seismrecord.dat',status='old')
         
		   do j=1,lsyn
                write(300,401) j*samp,(seis(i,j)+(i-1)*25.0,i=1,d)
		   enddo
401        format(1x,f6.1,18f10.3)
      









	end

!..........................................................

	subroutine  Rtcoef(vp,vp1,vp2,vp0,hs,r0,m,n,res) 
      implicit complex (a-h,o-z)
      complex res,r00,mm,nn
      real     vp,vp1,vp2,vp0,hs,r0,m,n 
           
           alpha3=cmplx(vp,0.) 
           alpha0=cmplx(vp0,0.)
           alpha1=cmplx(vp1,0.)
		   alpha2=cmplx(vp2,0.)      
           p=cmplx(hs,0.)
           r00=cmplx(r0,0.)
		   mm=cmplx(m,0.)
		   nn=cmplx(n,0.)

		   s11=alpha1*p
           s12=alpha2*p          
          
          c11=csqrt(1.0-s11**2)  !rushejiao  yuxuan
          c12=csqrt(1.0-s12**2)  !zheshejiao yuxuan
         
		  c13=csqrt((1.0-c11)/2.0) !入射角一半的yu玄
		  c14=csqrt((1.0+c12)/2.0) !折射角一半的yu玄
          c23=csqrt((1.0+c11)/2.0)
		  c24=csqrt((1.0-c12)/2.0)
		  ttas=c13*c14+c23*c24
         res=r0+mm*ttas**2+nn*(ttas/(csqrt(1-ttas**2))**2-ttas**2)
		return
        end


	 Subroutine wavelet(f,t0,w0,c,s,sampt,b,lb)
	 dimension  b(lb)
	  w   = w0/180.*3.1415927
      pi2 = 2.*3.1415927*f
	  a   = -f*f*alog(s)
      t   = -t0
      print*,'lb===',lb
      do  10  i = 1, lb
             b(i)= c * exp( a*t*t ) * cos( pi2*t+w )
             ccc=b(i)
10       t= t + sampt
       return
       end


	  Subroutine   fold(la,a,lb,b,lc,c,scaler)    
      dimension    a(la),b(lb),c(lc)                                
      call   zero(lc,c)
	  do   10   i=1,la
	  do   10   j=1,lb
	  k    = i+j-1
10	  c(k) = a(i)*b(j)+c(k)
      do   20   i=1,la+lb-1
20    c(i) = scaler*c(i)                                                                                              
      return                                                          
      end

	Subroutine   zero( l0, o )                                      
      dimension    o(l0)                                              
      do   10   i = 1, l0                                             
10    o(i)  = 0.                                                      
      return                                                          
      end