PROGRAM TWOBOX
	DIMENSION TITLE(4),PS1(51,30),PS2(51,30),PS12(51,30)
	COMMON N1X,N1Y,N2X,N2Y
c
c	************************************
c	Open a file to store the output
c	************************************        
c
	open(22,file='twobox.out')
	WRITE(6,1001)
c
c	******************************************************
c	The equetion variables : IOP(Wavefunction Number(1,2))
c	N1X,N1Y (quantum numbers for wavefunction 1)
c	******************************************************
c
	PRINT*,'GIVE ME IOP N1X N1Y'
	READ(5,*)IOP,N1X,N1Y
	PRINT*,IOP,N1X,N1Y
1001	FORMAT(1X,'PARTICLE IN A TWO DIMENSIONAL BOX',/,1X,107(1H*))
	PI=3.1416
	IF(IOP.EQ.2)READ(5,*)N2X,N2Y

c	**************************
c	Box dimentions X=51, Y=30
c	**************************

9	DO 2 I=1,51
	X=(I-1)/50.0
	DO 2 J=1,28
	Y=(J-1)/27.0
c
c	*********************
c	PS1 = Wavefunction 1
c	*********************
c
	PS1(I,J)=SIN(N1X*X*PI)*SIN(N1Y*Y*PI)*2.0
	IF(IOP.EQ.1)GO TO 3
c
c	********************
c	PS2 = Wavefunction 2
c	********************
c
	PS2(I,J)=SIN(N2X*X*PI)*SIN(N2Y*Y*PI)*2.0
c
c	******************************
c	PS12 = Product of PS1 and PS2
c	******************************
c
	PS12(I,J)=PS1(I,J)*PS2(I,J)
c
c	*********************************************
c 	Store the calculated data, in the output file
c	*********************************************
c
3       write(22,123)x,y,ps1(i,j),ps2(i,j),ps12(i,j)
2	CONTINUE
123     format(1x,5(f10.5,2x))
	CALL MAP(PS1,1)
1002	FORMAT(4I5)
	WRITE(6,103)
103	FORMAT(1H1)
	IF(IOP.EQ.1)GO TO 4
        WRITE(6,103)
	CALL MAP(PS2,2)
	WRITE(6,103)
	CALL MAP(PS12,3)
	CALL ING(PS12)
	WRITE(6,103)
	READ(5,1002)N1X,N1Y,N2X,N2Y
	IF(N1X.EQ.0)STOP
	GO TO 9
4       READ(5,1002)N1X,N1Y
	IF(N1X.EQ.0)STOP
	GO TO 9
	END

c	***********************************
c	Graphical representation of results
c	*********************************** 
	
	SUBROUTINE MAP(PS,IOUT)
	DIMENSION PS(51,30),IXY(51,30)
	COMMON NX1,NY1,NX2,NY2
	IF(IOUT.GT.1)GO TO 1
	WRITE(6,100)NX1,NY1
	GO TO 3
1	IF(IOUT.GT.2)GO TO 9
	WRITE(6,100)NX2,NY2  
	GO TO 3
9	CONTINUE
	WRITE(6,101)NX1,NY1,NX2,NY2
3	CONTINUE
100	FORMAT(1X,'MAP OF PSI(',I2,1X,I2,')*49.5',/)
101	FORMAT(1X,'MAP OF PSI(',I2,1X,I2,')*PSI(',I2,1X,I2,')*24.75',/)
	WRITE(6,105)
105	FORMAT(1X,86(1H*))
	IF(IOUT.EQ.3)GO TO 5
	DO 6 I=1,51
	DO 6 J=1,28
6	IXY(I,J)=NINT(PS(I,J)*49.5)
	GO TO 7
5	DO 8 I=1,51
	DO 8 J=1,28
8	IXY(I,J)=NINT(PS(I,J)*24.75)
7	CONTINUE
	DO 10 J=1,28
10	WRITE(6,107)(IXY(I,J),I=1,51)
107	FORMAT(1X,1H*,51I3,1H*)
	WRITE(6,105)
	RETURN
	END
	
	SUBROUTINE ING(PS)
	DIMENSION PS(51,30)
	A=1.0/51
	B=1.0/28
	DV=A*B
	SUM=0.0
	DO 1 I=1,51
	DO 1 J=1,28
c
c	*******************************
c	Calculation of the integral XX
c	*******************************
c
1	SUM=SUM+PS(I,J)
	XX=SUM*DV
	WRITE(6,100)XX
100	FORMAT(1X,'THE INTEGRAL OF THIS PRODUCT OVER ALL THE BOX SPACE=',
     1  1F10.3)          
	RETURN
	END