	PROGRAM PHSTAM
C	Program to calculate means and standard deviations of phase 
C	distributions of rgith & left ventricles and of two halves
C	of the left ventricle. The significance of the difference
C	in means in each case is calculated.
C	The program also sets up a histogram of the phase distribution
C	of the LV ,stores it in frame 3 of the multiple static phase
C	file, and calculates the FWHM and right and left skews of the
C	distribution.
C	The program provides an option to set up a playback buffer to
C	help visualise the sequence of contraction.
	VIRTUAL MATA(64,64),MAT(64,64),ROI(64,64)
	VIRTUAL CMAT(64,64)
	INTEGER*2 ROI
	INTEGER*2 SADMIN(256),P(75)
	INTEGER*2 CU1(64),CMAX,CU9,SKR,SKL
	LOGICAL*1 BMAT1(4096),BMAT(64,64)
	LOGICAL*1 NAME(14),L1,L2,LG(35),HEAD(29),M1(14),M2(14),M3(14)
	LOGICAL*1 RADMIN(510),ADMIN(512),CMAT,NAME1(14),LADMIN(512),ANS,BANS
	LOGICAL*1 TITLE(6)
	REAL*4 B(3,42),MEANRV,MEANLV,MEANL,MEANR,MEANLL,MEANLR
	EQUIVALENCE(BMAT,SADMIN,B,LADMIN),(BMAT1(2048),ADMIN(512))
	EQUIVALENCE (BMAT,BMAT1),(BMAT1(3000),RADMIN)
	DATA VIEW3/'HIST'/
	DATA CA/'CA  '/
	DATA HEAD/'P','H','A','S','E',' ','D','I','S','T','R','I','B',
	1'U','T','I','O','N',' ','P','A','R','A','M','E','T','E','R','S'/
	DATA M1/'F','W','H','M','=',' ',' ',' ',' ',' ',' ',' ',' ','%'/
	DATA M2/'R','I','G','H','T',' ','S','K','E','W','=',' ',' ','%'/
	DATA M3/'L','E','F','T',' ',' ','S','K','E','W','=',' ',' ','%'/
	DATA NAME1/'S','V','A',':','S','V','A','R','0','2','.','S','Y','S'/
	DATA TITLE/'P','H','S','T','A','M'/
	CALL CNTXSW
	CALL GAMSTP(SADMIN,TITLE)
	LP = 6
C
C	OPEN SVAR06 TO GET NUMBER OF FRAMES IN STUDY
C
	NAME1(10)= '6'
	CALL ASSIGN(1,NAME1,14,'OLD','NC')
	DEFINE FILE 1(53,256,U,N1)
	N1 = 1
	READ(1'N1)SADMIN
	CALL CLOSE(1)
	NFR =SADMIN(107)
C
C	phase image with  right and left ventricle ROIs
C	in SA2	
C
	NAME1(10) = '2'
	CALL ASSIGN(1,NAME1,14,'OLD','NC')
	DEFINE FILE 1(33,256,U,N1)
	N1=1
	READ(1'N1)SADMIN
C	read name of MS file from admin block of SA2
	CALL R50ASC(3,SADMIN(1),NAME(1))
	CALL R50ASC(3,SADMIN(2),NAME(5))
	CALL R50ASC(3,SADMIN(3),NAME(8))
	CALL R50ASC(3,SADMIN(4),NAME(12))
	NAME(4)=':'
	NAME(11)='.'
	CALL ASSIGN(2,NAME,14,'OLD',NC)
C	WRITE(7,299)(NAME(I),I=1,14)
299	FORMAT(' SECOND FILE',14A1)
	DEFINE FILE 2(52,256,U,N2)
	N1=2
C	read matrix data of phase image from SA
	DO 40 I=1,64,4
40	READ(1'N1)((MAT(K,J),J=1,64),K=I,I+3)
C	read ROI information from SA
	N1=18
	DO 50 I=1,64,4
50	READ(1'N1)((ROI(K,J),J=1,64),K=I,I+3)
C	read admin parameters
	N2=1
	READ(2'N2)ADMIN
	CALL FGADM1(ADMIN,P,B)
	WRITE(LP,604)
604	FORMAT('1')
	WRITE(LP,605)((B(I,J),I=1,3),J=1,2)
	WRITE(LP,605)(B(I,5),I=1,3)
	WRITE(LP,605)(B(I,9),I=1,3)
605	FORMAT(2X,6A4)
C	calculate parameters for RV
	SUMRV=0.0
	ITR=0
	DO 60 I=1,64
	DO 60 J=1,64
	IR=ROI(I,J)
	IF(IR.EQ.2.OR.IR.EQ.0)GO TO 60
	IF(MAT(I,J).EQ.0)GO TO 60
	SUMRV=MAT(I,J)+SUMRV
	ITR=ITR+1
60	CONTINUE
	RMNRV=SUMRV/ITR
	MEANR=(RMNRV/6283.)*360.
	MEANRV=RMNRV
	SUM=0.0
	DO 70 I=1,64
	DO 70 J=1,64
	IR=ROI(I,J)
	IF(IR.EQ.2.OR.IR.EQ.0)GO TO 70
	IF(MAT(I,J).EQ.0)GO TO 70
	DIFF=MAT(I,J)-MEANRV
	SUM=(DIFF*DIFF)+SUM
70	CONTINUE
	STDRV=SQRT(SUM/ITR)
	STDRV=(STDRV/6283.)*360.
	WRITE(LP,72)MEANR,STDRV
72	FORMAT(//2X,'MEAN RV,DEGREES=',F6.0,4X,'STD RV=',F6.0)
C	calculate parameters for LV
	SUMLV=0.0
	ITL=0
	DO 80 I=1,64
	DO 80 J=1,64
	IR=ROI(I,J)
	IF(IR.EQ.2.OR.IR.EQ.3)GO TO 82
C	zero matrix elements which are not in LV
	MAT(I,J)=0
	GO TO 80
82	IF(MAT(I,J).EQ.0)GO TO 80
	SUMLV=MAT(I,J)+SUMLV
	ITL=ITL+1
80	CONTINUE
	RMNLV=SUMLV/ITL
C	convert to degrees
	MEANL=(RMNLV/6283.)*360.
	MEANLV=RMNLV
	SUM=0.0
	DO 90 I=1,64
	DO 90 J=1,64
	IF(MAT(I,J).EQ.0)GO TO 90
	DIFF=MAT(I,J)-MEANLV
	SUM=(DIFF*DIFF)+SUM
90	CONTINUE
	STDLV=SQRT(SUM/ITL)
	STDLV=(STDLV/6283.)*360.
	WRITE(LP,92)MEANL,STDLV
92	FORMAT(/2X,'MEAN LV,DEGREES=',F6.0,4X,'STD LV=',F6.0)
500	DMEAN=MEANL-MEANR
	WRITE(LP,98)DMEAN
98	FORMAT(/2X,'DIFFERENCE IN MEANS=',F6.0)
C	CALL SIG(MEANR,MEANL,ITR,ITL,STDRV,STDLV)
C	calculate means,stdvns, and significance of difference in
C	means for left & right halves of LV
C	left boundary of LV
C	J1=LADMIN(369)/2
C	right boundary of LV
C	J2=LADMIN(371)/2
C	TYPE *,'J1',J1,'J2',J2
C	IX1=J1
C	IX2=(J2-J1)/2
C	CALL MEAN(IX1,IX2,MEANLL,STDLL,ITLL,MAT)
C	IX1=(J2-J1)/2
C	IX2=J2
C	CALL MEAN(IX1,IX2,MEANLR,STDLR,ITLR,MAT)
C	WRITE(LP,800)MEANLL,STDLL
800	FORMAT(//2X,'MEAN LEFT LV=',F6.0,'STD LEFT LV=',F6.0)
C	WRITE(LP,810)MEANLR,STDLR
810	FORMAT(/2X,'MEAN RIGHT LV=',F6.0,'STD RIGHT LV=',F6.0)
C	CALL SIG(MEANLL,MEANLR,ITLL,ITLR,STDLL,STDLR)
C set up histogram of LV distribution
510	SCALE=6283./64.
	DO 202 I=1,64
	CU1(I)=0
202	CONTINUE
	CMAX=0
	DO 305 I=1,64
	DO 305 J=1,64
	IF(MAT(I,J).EQ.0)GO TO 305
	IX=INT(FLOAT(MAT(I,J))/SCALE)
	IF(IX.EQ.0)GO TO 305
	CU9=CU1(IX)+1
	CU1(IX)=CU9
	IF(CMAX.GE.CU9)GO TO 305
	CMAX=CU9
	IMX=IX
305	CONTINUE
99	FORMAT(1X,F10.4)
	STEP=FLOAT(CMAX)/64.
	DO 332 I=1,64
	CU1(I)=FLOAT(CU1(I))/STEP
332	CONTINUE
	DO 360 I=1,64
	DO 360 J=1,64
360	MATA(I,J)=0
	DO 380 J=1,64
	IX=CU1(J)
	DO 380 I=1,IX
	MATA(I,J)=J
D	TYPE *,J,I,MAT(I,J)
380	CONTINUE
C calculate fwhm and skew from phase distribution
	DO 45 J=1,64
	IF(MATA(64,J).NE.0)GO TO 42
45	CONTINUE
42	KL=J
D	TYPE *,'KL',KL
	DO 33 J=KL,1,-1
	IF(MATA(32,J).NE.0)GO TO 33
	J1=J
	GO TO 34
33	CONTINUE
34	DO 37 L=KL,64
	IF(MATA(32,L).NE.0)GO TO 37
	J2=L
	GO TO 38
37	CONTINUE
38	IFWHM=((J2-J1)*100)/64
	SKR=((J2-IMX)*100)/64
	SKL=((IMX-J1)*100)/64
	WRITE(LP,606)HEAD
606	FORMAT(//2X,29A1)
	WRITE(LP,600)IFWHM
600	FORMAT(2X,'FWHM=',I3,'%')
	WRITE(LP,610)SKR
610	FORMAT(2X,'RIGHT SKEW=',I3,'%')
	WRITE(LP,620)SKL
620	FORMAT(2X,'LEFT SKEW=',I3,'%'//)
	DMEAN=ABS(DMEAN)
	WRITE(7,615)
615	FORMAT(//,'$IS THIS A RESTING STUDY- Y OR N?   ')
	READ(5,1000)ANS
1000	FORMAT(A1)
	IF(ANS.EQ.'Y')GO TO 1010
	DIS=STDLV+DMEAN*2.3+IFWHM*0.8
	GO TO 1020
1010	DIS=STDLV+DMEAN*3.3+IFWHM*4.6
	IDIS=INT(DIS)
	TYPE *,'  DISCRIMINANT FUNCTION VALUE',IDIS
	TYPE *,'  NORMAL<86'
	IF(DIS.GT.85.)GO TO 97
	GO TO 1030
1020	IDIS=INT(DIS)
	TYPE *,'  DISCRIMINANT FUNCTION VALUE',IDIS
	TYPE *,'  NORMAL<52'
	IF(DIS.GT.51.)GO TO 97
1030	TYPE *,'  NORMAL'
	GO TO 520
97	TYPE *,'  PROBABLY ABNORMAL'
C change admin parms for frame 3
520	N2=1
	READ(2'N2)ADMIN
	CALL FGADM1(ADMIN,P,B)
	B(1,10)=VIEW3
	P(1)=3
	P(2)=-34
	P(3)=52
	P(4)=-35
	P(5)=0
	P(6)=1
	P(7)=-17
	P(8)=17
	P(9)=3
	P(13)=1
	P(21)=3
	CALL FGADM2(ADMIN,P,B)
C	store frame 3 in multiple static file
	N2=36
	WRITE(2'N2)ADMIN
	DO 400 I=1,64,4
400	WRITE(2'N2)((MATA(K,J),J=1,64),K=I,I+3)
C convert histogram parameters to ascii for storing in comment block
	DO 401 I=1,34
401	LG(I)="000
	L1="01
	L2="200
	CALL ASC(IFWHM,M1)
	CALL ASC(SKR,M2)
	CALL ASC(SKL,M3)
	WRITE(2'2)L1,M1,L2,LG,L1,M2,L2,LG,L1,M3,L2
	CALL CLOSE(2)
902	WRITE(7,903)
903	FORMAT(//,'$DO YOU WANT PB BUFFER- Y OR N?   ')
	READ(5,904)L1
904	FORMAT(A1)
	IF(L1.EQ.'N')GO TO 302
	IF (L1.NE.'Y') GOTO 902
	CALL PBCRT(MAT,CMAT,BMAT,ADMIN,RADMIN,NAME,NFR,N1)
302	CALL CLOSE(1)
	REWIND LP
	CALL GAMCA(SADMIN,-1)
	END
