      PROGRAM WAVEPA
C
C     -----------------------------------------------
C     AUTHOR: ROBERT EHRLICH (GEORGE MASON U    1973)
C     Revised: Peter Signell (Mich. St. U.      3/86)
C     Revised: Gary Herzenstiel (Mich. St. U.   3/86)
C     Revised: E.J.D.Kales (Mich. St. U.       11/87)
C     Revised: Peter Signell (Mich. St. U.      3/02)
C
C     This program displays interface patterns which are
C     calculated for two point sources of waves.  The wave
C     amplitude from each source is assumed to be
C     independent of angle, and it varies as the inverse
C     first power of the distance from the source.
C
C     SUBROUTINES:
C          GETINF  - get values for DEL,RATIO,L1,L2,V from
C                    user
C          CALC    - Main program for calcualting the
C                    patterns
C          DOCOLS  - loops over 101 columns for a given row
C          DISLIN  - displays a line of the graphic output
C          QUERY   - asks if the user wants to do another
C
C                            W A V E P A T
C                                  :
C               --------------------------------------
C               :                  :                 :
C            GETINF               CALC             QUERY
C      (DEL,RATIO,L1,L2,V)  (DEL,RATIO,L1,L2,V)   (AGAIN)
C                                  :
C                        ------------------------
C                        :                      :
C                     DOCOLS                 DISLIN 
C          (DEL,RATIO,L1,L2,VS1,S2,J,RINT)    (RINT)
C
C     ***** DECLARATIONS *****
      REAL DEL,RATIO,L1,L2,V
      CHARACTER AGAIN
      OPEN(UNIT=3,FILE='(T:1041)m353p1f.out')
C                    
C     ***** M A I N   B O D Y *****
C
 10   CALL GETINF (   DEL,RATIO,L1,L2,V)
      CALL CALC   (DEL,RATIO,L1,L2,V   )
      CALL QUERY  (               AGAIN)  
      IF(AGAIN.EQ.'Y') THEN
          WRITE(3,'(A1)') 12
          GOTO 10
      ENDIF
      IF(AGAIN.EQ.'y') THEN
          WRITE(3,'(A1)') 12
          GOTO 10
      ENDIF
      STOP
      END
C
C     ***** G E T I N F *****
C
      SUBROUTINE GETINF(   DEL,RATIO,L1,L2,V)
C 
C     gets values for DEL,RATIO,L1,L2, & V from user.
C
C     ***** DECLARATIONS *****
      REAL DEL,RATIO,L1,L2,V
C
C     DEL   = separation between sources
C     RATIO = ratio of two source strengths
C     L1,L2 = wavelengths
C     V     = velocity of source 1. (waves V is assumed 1.0)
C
      PRINT *,'Please supply DEL, RATIO, L1, L2 & V values:'
      WRITE(*,20) 'DEL='
      READ *,DEL
      WRITE(*,20) 'RATIO='
      READ *,RATIO
      WRITE(*,20) ' L1='
      READ *,L1
      WRITE(*,20) ' L2='
      READ *,L2
      WRITE(*,20) '  V='
      READ *,V
  20  FORMAT(A4/)
      WRITE(*,24) DEL,RATIO,L1,L2,V
      WRITE(3,24) DEL,RATIO,L1,L2,V
  24  FORMAT('m353p1f.for    Wave Interference Patterns'/
     +'DEL=',F7.4,' R=',F7.4,' L1=',F7.4,' L2=',F7.4,/
     +' V=',F7.4)
      RETURN
      END
C
C     ***** C A L C *****
C
      SUBROUTINE CALC(DEL,RATIO,L1,L2,V   )
C       
C     shell for the calculation loops over all rows and
C     each column within a row, then prints the row graphic
C
C     ***** DECLARATIONS *****
      REAL DEL,RATIO,L1,L2,V,S1,S2,RINT(101)
      REAL XSCLE,YSCLE
      INTEGER J
C
      S1=1.0
      S2=RATIO
      XSCLE=0.78
      YSCLE=0.78
      XSCLE=INT(51.0*XSCLE)/51.0
      YSCLE=INT(31.0*YSCLE)/31.0
      WRITE(*,28) XSCLE,YSCLE
      WRITE(3,28) XSCLE,YSCLE
 28   FORMAT('XSCLE=',F7.5,'   YSCLE=',F7.5/)
C
C     loop over 53 ROWS and call subroutine to loop over
C     columns, then print the graphic line before doing the
C     next row:
C
      DO 30,J=1,53
           CALL DOCOLS(DEL,RATIO,L1,L2,V,S1,S2,J,XSCLE,YSCLE,RINT)
           CALL DISLIN(RINT   )
 30   CONTINUE
      RETURN
      END
C
C     ***** D O C O L S *****
C
      SUBROUTINE DOCOLS(DEL,RATIO,L1,L2,V,S1,S2,J,XSCLE,YSCLE,RINT)     
C
C     loops over 80 columns per row calculating combined waves intensity
C
C     ***** DECLARATIONS *****
      REAL DEL,RATIO,L1,L2,V,S1,S2,AJ,AK,X,Y,R1,R2,D
      REAL A,B,C,DISC,X0,A1,A2,PI,CONS,RINT(101)
      REAL XSCLE,YSCLE
      INTEGER K,J
C
      PI=3.14159265
CC      DO 40,K=1,101
      DO 40,K=1,80
C
C          (X,Y) are calculated from row and column indices
C          J and K
C
           AJ=J
           AK=K
           X=AK/(5.0*XSCLE)-51/5.0
           Y=AJ/(3.0*YSCLE)-31/3.0
C
C          R1 and R2 are the distances from each source to
C          (X,Y)
C
           R2=SQRT(X**2+(Y-DEL/2.0)**2)
           IF(V.EQ.0) THEN
C               the velocity of the source is ZERO:
                R1=SQRT(X**2+(Y+DEL/2.0)**2)
           ELSE
C               the velocity of the source is NON-ZERO:
                R2=L2
                IF(V.EQ.1.0) V=1.00001
                D=3.0*V
                IF(D.GT.10.0) D=10.0
C
C               solve using the quadratic formula:
C
                A=1.0-1.0/V**2
                B=2.0*D/V**2-2.0*X
                C=X**2+Y**2-(D/V)**2
                DISC=B**2-4.0*A*C
                R1=L1
                IF(DISC.GE.0) THEN
                   X0=(-B+SQRT(DISC))/(2.0*A)
                   IF(X0.LE.D) THEN 
                      R1=SQRT((X-X0)**2+Y**2)
                   ENDIF
                ENDIF
           ENDIF
C
C          add .000001 to R1, R2, and L2 to avoid division
C          by ZERO
C
           R1=R1+.000001
           R2=R2+.000001
           L2=L2+.000001
C
C          A1 & A2 are the amplitudes of each wave at point
C          (X,Y)
C
           A1=S1*SIN(2.0*PI*R1/L1)/R1
           A2=S2*SIN(2.0*PI*R2/L2)/R2
C
C          RINT is the resultant intensity of the two waves
C          combined, CONS is an arbitrary normalization
C          constant
C
           CONS=2.0/(S1**2+S2**2)
           RINT(K)=CONS*(A1+A2)**2
C
C          if the point (X,Y) is closer than 0.7 units to
C          either source, set the intensity digit to 0, to
C          make a circular blank around the location of each
C          source.
C
           IF(V.EQ.0)THEN
                IF(R1.LE.0.7) RINT(K)=0.0
                IF(S2.NE.0.AND.R2.LE.0.7) RINT(K)=0.0
           ENDIF
C
 40   CONTINUE
      RETURN
      END
C
C     ***** D I S L I N *****
C
      SUBROUTINE DISLIN(RINT   )
C
C     display a line of graphic output 
C
C     ***** DECLARATIONS *****
      REAL RINT(101),THRES(7)
      INTEGER L,M,FLAG,INDEX
      CHARACTER E*1041,PCHAR(7)
      DATA THRES/.005,.01,.02,.04,.08,.16,.32/
      DATA PCHAR/' ','I','-','X','O','#','M'/
C
C     each line of output is overprinted up to 7 times if
C     necessary
C
      INDEX=0
      DO 55,L=1,7
C
           DO 51,M=1,80
                INDEX=INDEX+1
C
C               if intensity exceeds threshold value
C               THRES(L), use character PCHAR(L), otherwise
C               use a blank.
C
                IF(RINT(M).GE.THRES(L)) THEN
                     E(INDEX:INDEX) = PCHAR(L)
                ELSE
                     E(INDEX:INDEX) = ' '
                ENDIF
 51        CONTINUE

           IF(L.EQ.7) GOTO 55
           
           DO 52 M=1,80
              INDEX=INDEX+1
              E(INDEX:INDEX) = CHAR(8)
 52        CONTINUE

 55   CONTINUE
      WRITE(3,54) E(1:1040)

 54   FORMAT(A1040)

      RETURN
      END
C
C     ***** Q U E R Y *****
C
      SUBROUTINE QUERY(    AGAIN)
C 
C     ***** DECLARATIONS *****
      CHARACTER AGAIN
C
      PRINT *
      WRITE(*,62)
 62   FORMAT(
     +'Would you like to do this again using new values?'/
     +'(Y/N)? ')
      READ( *,'(A1)')AGAIN
      RETURN
      END
C     ***** P R O G R A M    E N D *****
