Ask Question

Name:
Title:
Your Question:

Answer Question

Name:
Your Answer:
User Submitted Source Code!


Description:
  compile
Language: FORTRAN
Code:
C     Last change:  OR   21 Jun 2005    1:29 pm

C     SPECIFY THE PROBLEM SPACE SIZE HERE.  NX, NY AND NZ SPECIFY
C     THE SIZE OF THE PROBLEM SPACE IN THE X, Y AND Z DIRECTIONS.

c      PARAMETER (NPMLS=5,NXS=109,NX1S=NXS-1,ISRCS=55)
      PARAMETER (NPMLS=9,NXS=117,NX1S=NXS-1,ISRCS=59)
c      PARAMETER (NPMLS=11,NXS=121,NX1S=NXS-1,ISRCS=61)
c      PARAMETER (NPMLS=17,NXS=133,NX1S=NXS-1,ISRCS=67)
c      PARAMETER (NPMLS=33,NXS=165,NX1S=NXS-1,ISRCS=83)
c      PARAMETER (NPMLS=33,NXS=1065,NX1S=NXS-1,ISRCS=533)

c      PARAMETER (NPMLS=9,NXS=617,NX1S=NXS-1,ISRCS=309)
c      PARAMETER (NPMLS=33,NXS=1065,NX1S=NXS-1,ISRCS=533)

c      PARAMETER (NPMLS=5,NXS=7409,NX1S=NXS-1,ISRCS=405)
c      PARAMETER (NPMLS=9,NXS=7417,NX1S=NXS-1,ISRCS=409)
c      PARAMETER (NPMLS=17,NXS=7433,NX1S=NXS-1,ISRCS=417)
c      PARAMETER (NPMLS=33,NXS=7465,NX1S=NXS-1,ISRCS=433)
c      PARAMETER (NPMLS=65,NXS=7529,NX1S=NXS-1,ISRCS=465)
c       PARAMETER (NPMLS=65,NXS=17529,NX1S=NXS-1,ISRCS=10465)


      PARAMETER (NSTOP=01.*16*1024)

c      PARAMETER (NSTOP=4)

      PARAMETER (eps0=8.854e-12, xmu0=12.566e-7, eta0=376.734)

c      PARAMETER (cond=0.620,epsinf=50.,epss=160,t0=5.88e-9)
c      PARAMETER (cond=0.0,epsinf=1.,epss=2.25,f0=.6332e16,delta0=.28e16
c     1 ,alphan=1,x03=0)

c      PARAMETER (cond=0.0,epsinf=2.25,epss=5.25,fl=63.7e12,delta0l=1.e9
c     1,fr=1.396e13,delta0r=3.1245e13
c     1,alphan=0.7*1,x03=0.07*1)

c      PARAMETER (cond=0.0,epsinf=1.,epss=2.25,fl=.6332e16,delta0l=.28e16
c     1,fr=1.396e13,delta0r=3.1245e13
c     1,alphan=0.7*0,x03=0.07*0)

c      PARAMETER (cond=0.0,epsr=1.,epss=3.,f0=10.e9,delta0=3.e8)

c      PARAMETER (epsinf=1.,epss=2.,cond=0)
c      PARAMETER (fp=1.1027e9,delta0=0,f0=0.15915e9)


      PARAMETER (cond=0.,epsr=1.)

c      PARAMETER (cond=0.0,epsinf=1.,epss=2.25,fl=.6332e16,delta0l=.28e16
c     1,fr=1.396e13,delta0r=3.1245e13
c     1,alphan=0.7*0,x03=0.07*0)


C     DEFINE NSTOP, (MAXIMUM NUMBER OF TIME STEPS) HERE.
C
      COMMON/EfldS/ezs(nxs),ez1s(nxs)
     1             ,fzxs(nxs),fzx1s(nxs)
     1            ,dzs(nxs),dz1s(nxs)
     1             ,pzxs(nxs),pzx1s(nxs)
     1            ,szs(nxs),sz1s(nxs)
     1            ,izs(nxs),iz1s(nxs)
     1            ,sz11s(nxs),sz21s(nxs)
     1            ,szls(nxs),szl1s(nxs)
     1            ,szl11s(nxs),szl21s(nxs)
     1            ,szrs(nxs),szr1s(nxs)
     1            ,szr11s(nxs),szr21s(nxs)


      COMMON/HfldS/hys(nxs),hy1s(nxs)
     1            ,gyxs(nxs),gyx1s(nxs)
     1            ,bys(nxs),by1s(nxs)
     1            ,sys(nxs),sy1s(nxs)
     1            ,iys(nxs),iy1s(nxs)
     1            ,sy11s(nxs),sy21s(nxs)


      COMMON/cons_tri/x(nxs),y(nxs),mgamma(nxs),g(nxs)
     1 ,aa1(nxs),bb1(nxs),cc1(nxs),d(nxs),x1(nxs),x2(nxs)
     1, aa2(nxs),bb2(nxs),cc2(nxs),perm(nxs)
C
      COMMON/EXTRAS/N,DT,T,NPTS,C,PI,DELX,DELY,EZ11,EZ21,dtdx,TT
     1,alpha0,alpha1,alpha2,alpha3,beta1,beta2,beta3,ez1,ez2,ez3
     1,ez4,ez5,cdtdx,dx,xx,yy,zz,ww,cons0,cons1,cons2,beta,gama,alpha
     1,ty,af,b0,b1,a1,t1,alphal,alphar,betal,gamal,betar,gamar,eps1

      COMMON/CONPML/C1(nPMls),C2(nPMls), C3(nPMls), C4(nPMls)
     &,sig1(npmls),sig2(npmls)
     1,C11(nPMLs), C12(nPMLs), C13(nPMLs), C14(nPMLs)
       
      COMMON/PMLCEFS/fi3s(nxs),gi3s(nxs),fi2s(nxs),gi2s(nxs)
     1,fi1s(nxs),gi1s(nxs)
     &,b0es(nxs),b1es(nxs),a1es(nxs),be0s(nxs),be1s(nxs)
     &,b0hs(nxs),b1hs(nxs),a1hs(nxs),bh0s(nxs),bh1s(nxs)


     C     Last change:  OR   22 Jun 2005    9:25 am

      INCLUDE 'adi1dc.for'
C
      OPEN(UNIT=10,FILE='adir1.dat',STATUS='UNKNOWN')
      OPEN(UNIT=20,FILE='afadi2.dat',STATUS='UNKNOWN')
      OPEN(UNIT=30,FILE='afadi3.dat',STATUS='UNKNOWN')
      OPEN(UNIT=40,FILE='afadi4.dat',STATUS='UNKNOWN')

c     ZERO PARAMETERS, GENERATE PROBLEM SPACE, INTERACTION OBJECT,
C     AND EXCITATION

      CALL ZEROS
c      CALL ZEROL
      CALL SETUP
C
C     *****************************************************************
C
C     MAIN LOOP FOR FIELD COMPUTATIONS AND DATA SAVING
C
C     *****************************************************************
      T=0.0
      ty=0
      
      DO 100 N=1,NSTOP
C
      WRITE (*,*) N
C
C     ADVANCE SCATTERED ELECTRIC FIELD
C
      CALL EZSFS1
      CALL HYSFS1


      CALL EZSFS2
      CALL HYSFS2

C     SAMPLE FIELDS IN SPACE AND WRITE TO DISK
C
      CALL DATSAV
C
       Ty=Ty+dt
 100  CONTINUE        

C     CLOSE DATA FILES
C
      CLOSE (UNIT=40)
      CLOSE (UNIT=30)
      CLOSE (UNIT=20)
      CLOSE (UNIT=10)
      STOP
      END
C     ********************************************************
      SUBROUTINE SETUP
C
      INCLUDE 'adi1dc.for'
C
C     THIS SUBROUTINE INITIALIZES THE COMPUTATIONS
C
      DT=25E-12

C     DEFINE PI DELX, AND C
C

        epsinf=epsr
      C=1.0/SQRT(XMU0*EPS0)
      PI=4.0*ATAN(1.0)
      DELX=1.*2.*C*DT
      DELY=DELX 
      dtdx=c/delx
      cdtdx=.5*c*dt/delx
      delx=1.e-3
      DELY=DELX
      dt=1*delx/(1.*c/SQRT(epsinf))
      dtdx=c/delx
      cdtdx=.5*c*dt/delx

C     CALCULATE DT--THE MAXIMUM TIME STEP ALLOWED BY THE
C     COURANT STABILITY CONDITION
C
      DTXI=C/DELX
      DTYI=C/DELY
c      DT=1./SQRT(DTXI**2+DTYI**2)

c        alpha=(1-dt/4/t0)/(1+dt/4/t0)
c        beta=(epsinf+epss*dt/4/t0)/(1+dt/4/t0)
c        gama=(epsinf-epss*dt/4/t0)/(epsinf+epss*dt/4/t0)

        ww=cond*dt/(4.*eps0)

c        cons0=(beta*gama-ww)/(beta+ww)
c        cons1=(alpha-1)/(beta+ww)
c        cons2=1./(beta+ww)

c        cons0=1
c        cons1=0.
c        cons2=1

        cons0=1
        cons1=0.
        cons2=1

      ww=cond*dt/(4.*eps0)

      cons0=(EPSR-ww)/(EPSR+ww)
      cons2=1./(EPSR+ww)
      epsinf=epsr


C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        nexp=4
     rmax=1.e-5
     xnpml=npmls-1
     c0=3e8
        npml=npmls
     do i=1,npml
       sigmax=-1*(nexp+1)*eps0*c0*log(rmax)/(2.*delx*xnpml*SQRT(epsinf))
        c1(i)=(.25*dt/eps0)*sigmax*((i-1.)/xnpml)**nexp
        c2(i)=(.25*dt/eps0)*sigmax*((i-.5)/xnpml)**nexp
        c3(npml+1-i)=c2(i)
        c4(npml+1-i)=c1(i)
     end do

     do idum=1,npml
        fi3s(idum)=(1.-c4(idum))/(1.+c4(idum))
        fi2s(idum)=2*c4(idum)*1./(1+c4(idum))
        gi3s(idum)=(1.-c3(idum))/(1.+c3(idum))
        gi2s(idum)=2*c3(idum)*1./(1.+c3(idum))
     end do
     icount=0
     do idum=nxs-npml,nxs
        icount=icount+1
        fi3s(idum)=(1.-c2(icount))/(1.+c2(icount))
        fi2s(idum)=2*c2(icount)*1./(1+c2(icount))
        gi3s(idum)=(1.-c1(icount))/(1.+c1(icount))
        gi2s(idum)=2*c1(icount)*1./(1.+c1(icount))
     end do

      RETURN
      END
C     ******************************************************************
      SUBROUTINE DATSAV
C
      INCLUDE 'adi1dc.for'

c      SUMM=0.0
c      DO 25 J=JSRCL-25+1,JSRCL+25-1
c      DO 25 I=ISRCL-50+1,ISRCL+50-1
c      C011=eel(I,J)
c      SUMM=SUMM+(C011-ees(I-150,J-175))**2
c25    CONTINUE
c      WRITE(*,*)SUMM
c      WRITE(10,*)SUMM

c      IF(N.EQ.100)THEN
c      DO 13 J=JSRCL-25+1,JSRCL-25+1
c      DO 13 I=ISRCL-50+1,ISRCL+50-1
c      C011=ezl(I,J)
c      C022=(-C011+ezs(I-150,J-175))
c      WRITE(20,*)C022
c13    CONTINUE
c      ENDIF

      WRITE(10,*)Ty,ezs(isrcs),ezs(isrcs-20),ezs(isrcs-49)
      WRITE(*,*)ezs(isrcs-49)


10    RETURN
      END
C     ******************************************************************
      FUNCTION EZIS(I)
C
      INCLUDE 'adi1dc.for'

C     THIS FUNCTION COMPUTES THE Z COMPONENT OF THE INCIDENT
C     ELECTRIC FIELD
C
      XN=N
      IF((XN*DT).LE.1.E-9)THEN
        ezis=(1.e-9/320.0)*(0.0*10.0
     &  -15.0*(-2.0*pi*1.0/1.e-9)*
     &sin(2.0*pi*1.0*xn*2.5e-11/1.e-9)
     &   +6.0*(-2.0*pi*2.0/1.e-9)*
     &sin(2.0*pi*2.0*xn*2.5e-11/1.e-9)
     &   -1.0*(-2.0*pi*3.0/1.e-9)*
     &sin(2.0*pi*3.0*xn*2.5e-11/1.e-9))

      EZIS=(1./320.)*(10.-15.*COS(XN*DT*2.*PI/1.E-9)
     1+6.*COS(XN*DT*2.*PI*2./1.E-9)
     1-COS(XN*DT*2.*PI*3./1.E-9))
      ELSE
      EZIS=0.0
      ENDIF

      t=xn*dt
      t00=210.*dt
      Tx=80.*dt

c      ezis=1000*EXP(-(t-t00)**2/Tx**2)


      fmax=1.e9
      fmin=0.0
      pw=2.*(sqrt(6.)/(pi*(fmax-fmin)))
      omeg=pi*(fmax+fmin)
      t00=4.*pw
      t=xn*dt

      EzIS=EXP(-((t-t00)/pw)**2)
c     $*SIN(2.*PI*FR*(XN*DT-42.*DT))
     $*SIN(omeg*t)

      WRITE(40,*)T,ezis


      RETURN
      END
C     ********************************************************
      SUBROUTINE EZSFS1
C
      INCLUDE 'adi1dc.for'
C
C     THIS SUBROUTINE UPDATES THE EZ SCATTERED FIELD COMPONENTS
C
      do i=2,nx1s
      d(i)=cons0*ezs(i)-cons1*dzs(i)+cons2*
     1(-1*gi3s(i)*fzxs(i)
     1+cdtdx*(1-gi2s(i))*(
     1 hys(i)-hys(i-1)
     1-(fi3s(i)*gyxs(i)-fi3s(i-1)*gyxs(i-1)))
     1)
      end do

      do i=2,nx1s
      aa2(i)=-cdtdx*(1-gi2s(i))*cdtdx*(1-fi2s(i-1))*cons2
      bb2(i)=1+cdtdx*cdtdx*(1-gi2s(i))*(2-fi2s(i)-fi2s(i-1))*cons2
      cc2(i)=-cdtdx*(1-gi2s(i))*cdtdx*(1-fi2s(i))*cons2
      end do

      i=2

      x2(i)=cc2(i)/bb2(i)
      g(i) = d(i)/bb2(i)

      do i=3,nx1s-1
      IF(i.eq.isrcs)then
     x2(i) =ezis(isrcs)
      else
        x2(i)=cc2(i)/(bb2(i)-aa2(i)*x2(i-1))
        g(i) =(d(i)-aa2(i)*g(i-1))/(bb2(i)-aa2(i)*x2(i-1))
      endif
      end do

      i=nx1s
      x2(i)=(d(i)-aa2(i)*g(i-1))/(bb2(i)-aa2(i)*x2(i-1))

      do i=nx1s-1,2,-1
      IF(i.eq.isrcs)then
      x2(i) =ezis(isrcs)
      else
      x2(i) = g(i) - x2(i) *x2(i+1)
      endif
      end do

      do i=2,nx1s
      IF(i.eq.isrcs)then
      ez1s(i) =ezis(isrcs)
      else
      ez1s(i)=x2(i)
      endif
      end do

      ez1s(isrcs) =ezis(isrcs)

      do i=2,nx1s
      dz1s(i)=alpha*dzs(i)+beta*(ez1s(i)-gama*ezs(i))
      end do

      do i=2,nx1s
      fzx1s(i)=gi3s(i)*fzxs(i)+cdtdx*gi2s(i)*(
     1hys(i)-hys(i-1)
     1-(fi3s(i)*gyxs(i)-fi3s(i-1)*gyxs(i-1))
     1+cdtdx*(1-fi2s(i))*ez1s(i+1)
     1-cdtdx*(2-fi2s(i)-fi2s(i-1))*ez1s(i)
     1+cdtdx*(1-fi2s(i-1))*ez1s(i-1))
      end do

C      ez1s(isrcs) =ezis(isrcs)

      RETURN
      END
C     ********************************************************
      SUBROUTINE HYSFS1
C
      INCLUDE 'adi1dc.for'
C
C     THIS SUBROUTINE UPDATES THE HY SCATTERED FIELD COMPONENTS
C
      DO 10 I=1,NX1S
      gyx1s(i)=fi3s(i)*gyxs(i)+cdtdx*fi2s(i)*(ez1s(i+1)-ez1s(i))
      hy1s(i)=hys(i)+cdtdx*(ez1s(i+1)-ez1s(i))-gyx1s(i)
 10   CONTINUE
C
      RETURN
      END
C     ********************************************************
C     ********************************************************
      SUBROUTINE EZSFS2
C
      INCLUDE 'adi1dc.for'
C
C     THIS SUBROUTINE UPDATES THE EZ SCATTERED FIELD COMPONENTS
C
      do i=2,nx1s
      ezs(i)=cons0*ez1s(i)-cons1*dz1s(i)+cons2*
     1(cdtdx*(hy1s(i)-hy1s(i-1))
     1-fzx1s(i)
     1)
      end do

      ezs(isrcs) =ezis(isrcs)

      do i=2,nx1s
      dzs(i)=alpha*dz1s(i)+beta*(ezs(i)-gama*ez1s(i))
      end do

      do i=2,nx1s
      fzxs(i)=gi3s(i)*fzx1s(i)+cdtdx*gi2s(i)*(
     1 hy1s(i)-hy1s(i-1))
      end do

      RETURN
      END
C     ********************************************************
      SUBROUTINE HYSFS2
C
      INCLUDE 'adi1dc.for'
C
C     THIS SUBROUTINE UPDATES THE HY SCATTERED FIELD COMPONENTS
C
      DO 10 I=1,NX1S
      gyxs(i)=fi3s(i)*gyx1s(i)+cdtdx*fi2s(i)*(ez1s(i+1)-ez1s(i))
      hys(i)=hy1s(i)+cdtdx*(ez1s(i+1)-ez1s(i))-gyx1s(i)
 10   CONTINUE
C
      RETURN
      END
C     ********************************************************
      SUBROUTINE ZEROS
C
      INCLUDE 'adi1dc.for'
C
C     THIS SUBROUTINE INITIALIZES VARIOUS ARRAYS AND CONSTANTS TO ZERO.
C
      T=0.0

      DO 10 I=1,NXS

      dzs(i)=0.0
      dz1s(i)=0.0

      fzxs(i)=0.0
      ezs(i)=0.0

      fzx1s(i)=0.0
      ez1s(i)=0.0

      gyxs(i)=0.0
      hys(i)=0.0

      gyx1s(i)=0.0
      hy1s(i)=0.0

      fi3s(i)=1.0
      fi2s(i)=.0
      gi3s(i)=1.0
      gi2s(i)=.0

 10   CONTINUE

      RETURN
      END
          
Comments: