C	PROGRAM DISASTER.FT
C	-------------------
C
C
C
	SUBROUTINE DISASTER(IH,FDELTA,JLENGTH,IDUMP)
C
C
C	P LEMKIN
C	NATIONAL CANCER INSTITUTE, DBD, IPU
C	NATIONAL INSTITUTES OF HEALTH
C	BETHESDA, MD
C	
C
C	JUNE 14, 1977
C	MARCH 15, 1977
C	MARCH 9, 1977
C	MARCH 8, 1977
C	JAN 11, 1977
C
C	PURPOSE
C	-------
C	COMPUTE THE DISASTER ANALYSIS CORRECTIONS FOR HISTOGRAM
C	IH[1:256] WHERE IH IS A DOUBLE PRECISION (EAE) ARRAY
C	DIMENSION IH(512). THE DISASTER ANALYSIS IS TAKEN FROM
C	A PAPER BY SMITH, SOLMON, AND WAGNER "PRACTICAL AND
C	MATHEMATICAL ASPECTS OF THE PROBLEM OF RECONSTRUCTING OBJECTS FROM
C	RADIOGRAPHS", GIVEN AT THE FAR WEST SECTIONAL MEETING OF THE
C	AMERICAN MATH. SOC., MONTEREY, CALIF. APRIL 1975.
C
C	
S	OPDEF DISP1 6435
S	OPDEF DISP2 6436
C
	DIMENSION IVAL(2)
	DIMENSION FJUNK(20),H(286),HNEW(256)
C
C	[1] CVT HIST TO F.P.
C	DEFINE DPCVRT DEFS
	LFTOD=+1
	LDTOF=-1
C
	KIHPTR=0
S	TAD \IH# /SAVE THE PTR
S	DCA \KIHPTR
C
	SUM=0.0
	DO 100 I=1,256
S	CLA CMA
S	TAD \I
S	DISP1
S	CLA
	IVAL=IH
S	INC \IH#
S	TAD I \IH
S	DCA \IVAL#
S	INC \IH#
C
	CALL DPCVRT(IVAL,A,LDTOF)
	H(I)=A
	SUM=SUM+A
100	CONTINUE
C
C	RESTORE THE PTR TO IH
S	TAD \KIHPTR
S	DCA \IH#
C
C
C	DEFINE THE SPOOLING SWITCH
S	TAD I \IDUMP
S	AND (1
S	DCA \IOUTSPOOL
C
C
C	[2] DEFINE THE WINDOW SIZE
	L=JLENGTH
C	EXPRESS DELTA AS %
	DELTA=0.001*FDELTA
C
C
C	[3] FIX UP THE END POINTS  "FOR NOW"
	DO 300 I=1,20
	FJUNK(I)=H(1)
300	H(256+I)=H(256)
C
C
C	[4] PERFORM THE DISASTER ANALYSIS
S	JMS CHECK
C
C
C	[5] CVT H BACK TO REFERENCE IH D.P. ARRAY
	DO 500 I=1,256
	CALL DPCVRT(IVAL,HNEW(I),LFTOD)
	IH=IVAL
S	INC \IH#
S	TAD \IVAL#
S	DCA I \IH
S	INC \IH#
500	CONTINUE
	RETURN
C
C
C
C	***********************************************************
C	*SUBROUTINE     C H E C K
C	***********************************************************
C	DISASTER CORRECTION
C
S	CPAGE 3
S RCHECK,	JMP I CHECK
S CHECK, 	0 /ENTRY
C
	DO 1010 IX=1,256
S	CLA CMA
S	TAD \IX
S	DISP2
S	CLA
C
C	COMPUTE: A
S	JMS A
C
C	EXPRESS FNOISE AS %
	FNOISE=100.0*(ABS(H(IX)-A)/SUM)
C
C
C	IF IDUMP AND FNOISE > 0
C		THEN LIST IT;
	IF(IDUMP)600,610,600
600	IF(FNOISE)610,610,611
611	LSNUM=1+IOUTSPOOL
	LSNEW=1
	DO 1611 INDEX=1,LSNUM
	WRITE(LSNEW,612)IX,H(IX),FNOISE,DELTA,ABS(FNOISE-DELTA)
612	FORMAT(' H(',I3,')=',F11.0,
     1', %NOISE(X,', DELTA=',F8.3,', /N(X)-DELTA/=',F8.3)
1611	LSNEW=3
C
C
C	IF |FNOISE| > DELTA
C		THEN HNEW(IX) <== APRIME(IX)
C		ELSE HNEW(IX)<==H(IX);
C
610	IF(FNOISE-DELTA)1012,1012,1011
C
C		THEN REPLACE HNEW(IX)
1011	CONTINUE
S	JMS APRIME
	HNEW(IX)=APRIME
	GOTO 1010
C
C	ELSE
1012	HNEW(IX)=H(IX)
1010	CONTINUE
S	JMP RCHECK
C
C	***********************************************************
C	*SUBROUTINE     A
C	***********************************************************
C	COMPUTE:
C		A=(1/L)*SUM[I=(X-L/2):(X+L/2) OF H(I)]
C
S	CPAGE 3
S RA,	JMP I A
S A, 	0 /ENTRY
C
	A=0.0
	L2=L/2
	FL=L
C
	LL=IX-L2
	LU=IX+L2
	DO 1000 I=LL,LU
1000	A=A+H(I)
C
	A=A/FL
S	JMP RA
C
C	***********************************************************
C	*SUBROUTINE     APRIME
C	***********************************************************
C	COMPUTE:
C		APRIME=(1/L)*SUM[I=(X-L/2)(X-1), (X+1):(X+L/2) OF H(I)]
C
S	CPAGE 3
S RAPRIME,	JMP I APRIME
S APRIME, 	0 /ENTRY
C
	APRIME=0.0
	L2=L/2
	FLM1=L-1
C
	LL1=IX-L2
	LU1=IX-1
	LL2=IX+1
	LU2=IX+L2
	DO 1001 I=LL1,LU1
1001	APRIME=APRIME+H(I)
C
	DO 1002 I=LL2,LU2
1002	APRIME=APRIME+H(I)
C
	APRIME=APRIME/FLM1
S	JMP RAPRIME
	END
    