$STORAGE: 2
C COPYRIGHT (C) 1983 GLENN EVERHART
c all rights reserved
C G H A S P   - Generalized Histogram And Scatter Plot
C REQUIREMENTS:
C COMMONS /EXTRA/ AND /PLOTS/ MUST EXIST, AND ARRAY MA'S DIMENSION (AS 4
C BYTE INTEGERS) MUST BE PLACED INTO NDLTY.
C NPLTS IS THE NUMBER OF PLOTS TO BE GENERATED; THEY ARE ALLOCATED OUT OF
C ARRAY MA DYNAMICALLY.
C THE VARIABLES IN THE /PLOTS/ COMMON HAVE THE FOLLOWING MEANINGS:
C NDIM IS THE NUMBER OF DIMENSIONS. 1 IS HISTOGRAM, 2 IS SCATTER PLOT
C XMIN,YMIN ARE X,Y MIN COORDS IN THE HISTOGRAM
C DX,DY ARE BIN SIZES
C NBINX,NBINY ARE NUMBER OF BINS IN X AND Y (NOTE GHASP WILL INDICATE NUMBER
C	OF OVERFLOWS)
C TITLE IS AN ARRAY OF CHARACTERS USED TO PRINT OUT THE TITLE FOR THE PLOT.
C THE SUBROUTINE INTERFACE IS TO CALL THE SUBROUTINE PLOT.
C
C CALL:
C   CALL PLOT(XVAL,YVAL,IFUNCT,NPLT)
C  WHERE XVAL AND YVAL ARE X,Y COORDINATES FOR THE PLOT IF SCATTERPLOT, OR
C X IS THE COORDINATE AND Y THE WEIGHT IF A HISTOGRAM.
C IFUNCT IS -1, 0, 1, OR 2.
C	-1 MEANS INITIALIZE; CALL PLOT ONCE THIS WAY TO SET UP THE NUMBER
C		OF HISTOGRAMS AND INITIALIZE ITS SCRATCH VARIABLES.
C	0 MEANS INITIALIZE VARIABLES FOR A GIVEN PLOT NUMBER. THIS EXPECTS
C		YOU HAVE SET THE PLOTS COMMON VARIABLES UP BEFORE THE
C		CALL. SET XVAL TO 4H/ DIM/ AT THIS CALL TO PRINT SOME
C		INFORMATION ABOUT HOW MUCH OF THE PLOT ARRAY IS USED UP;
C		THIS WILL ALLOW YOU TO CHANGE THE SIZE OF MA TO WHAT
C		IS REALLY NEEDED.
C	1 MEANS ENTER A POINT IN THE HISTOGRAM/SCATTER PLOT, USING THE X AND
C		Y VALUES. NOTICE THAT THE COMMON /PLOTS/ VARIABLES ARE NOT
C		NECESSARILY THE SAME AS AT IFUNCT=0 TIME; ONCE THE PLOT IS
C		INITIALIZED YOU JUST ADD POINTS AND PLOT.
C	2 MEANS PLOT THE HISTOGRAM OR SCATTER PLOT. THE XVAL ARGUMENT IS
C		IMPORTANT AT THIS TIME; PLOT NUMBER MUST BE GIVEN.
C
C A VARIETY OF OPTIONS FOR PLOT FORMAT EXIST AND ARE ENCODED BY THE LETTER
C USED IN THE XVAL ARGUMENT OF PLOT AT THE TIME YOU CALL IT WITH THE IFUNCT
C ARGUMENT OF 2. TWO OF THESE ARE THAT THE PLOT CAN BE MADE AS HIGH AS IT
C NEEDS TO BE TO PLOT THE DATA. THIS IS THE VARY COMMAND AND IS ENCODED AS
C 4H/   V/ (SEE EXAMPLE CALLER PROGRAM). ANOTHER OPTION IS THE HACK OPTION,
C CUTTING OFF THE PLOT AT ONE PAGE. THIS USES THE VALUE 4H/   H/. ONE CAN
C ALSO SCALE THE PLOT TO FIT ON A PAGE; 4H/   S/ WILL DO
C THIS. THERE ARE SOME DENSITY PLOTS AVAILABLE ALSO FOR SCATTER PLOTS;
C THE NORMAL PLOTS ARE 2 DIGIT NUMBERS (6 BITS ARE USED, PACKED 5 BINS TO
C A WORD, FOR COUNTING NUMBERS PER BIN).
C USE 4H/   Q/ FOR SHADED SCATTER PLOTS; DENSITY WILL BE APPROXIMATE ONLY
C BUT GHASP WILL ATTEMPT TO PLOT A SCATTER PLOT SHADED. NOTE A PLOT CAN
C BE PRINTED OUT MORE THAN ONCE, IN DIFFERENT FORMATS, SO A PLOT MAY BE
C PLOTTED NUMERICALLY WITH 4H/   V/ AND IN SHADED MODE AS WELL.
C
C COPYRIGHT (C) 1983, 1984 GLENN AND MARY EVERHART
C ALL RIGHTS RESERVED
C	PCC GRAF
C  NOTE: REQUIRES LUN 4 FOR TERMINAL OUTPUT AND LUN 5 FOR
C  TERMINAL INPUT.
C (GHASP USES 6 FOR LP:)
C
C GRAPHICS INTERFACE AND OUTPUT FROM PCC SPREADSHEET
C
C  GLENN EVERHART, 23-JAN-83
C
C SYNTAX AND USAGE:
C
C  This program is designed to allow an interactive user to enter
C a single command line to the program which it will parse (using
C the special version of VARSCN in GVARSCN) and allow graphic output
C from PortaCalc saved spread sheets. The assumption made is that
C the sheet has been saved with the PPN or PDN command. The filename
C must appear in the command line and variables in the file (named as
C though the cursor had been in cell A1 when the PPN / PDN was done)
C may be histogrammed or scatterplotted against each other.
C
C The GHASP routine (a FORTRAN plot package for ordinary printers)
C will be used for this version.
C
C Input syntax:
C  NN or LL filename.ext V1:V2 c V3:V4 +switches
C
C  where
C
C  an L in columns 1 or 2 takes log of 1st or 2nd range numbers (base 10),
C
C  filename appears at the start of the command line after a space
C  and with a space following it and is a valid RSX or VMS file spec.
C
C  V1:V2 and V3:V4 are ranges. V3:V4 is optional and its presence implies
C  a scatter plot. These ranges must be either a row or a column or part
C  of them. If only range V1:V2 is present, a histogram will be done using
C  the Scale option of GHASP to fit the plot onto a page. The plot will be
C  set up for 100 bins horizontal, 50 vertical.
C    If the V3:V4 range exists, the character Q in the c position (the
C  are required) will result in a "density" plot in which the program will
C  attempt to print darker in filled bins. This is crude and the default is
C  to use a 2 digit number. Again, plot size will be scaled to 50 by 50
C  bins.
C
	CHARACTER*1 LINE(128),KLET,LLET
	INTEGER*4 NDLTY,NPLTS,MA(1000)
	INTEGER*4 NDIM,TITLE(19)
	CHARACTER*1 LTITL(76)
	CHARACTER*1 LLA,LLB
	EQUIVALENCE(TITLE(1),LTITL(1))
	INTEGER*2 IXTR
	COMMON/IXTR/IXTR
	INTEGER*4 NBINX,NBINY
	REAL*4 XMIN,YMIN,DX,DY
	INTEGER*2 KK,LS1,LS2,LQ
	REAL*4 VEC1(300),VEC2(300)
	CHARACTER*1 IONM(50)
	COMMON/EXTRA/NDLTY,NPLTS,MA
	EXTERNAL INDEX
	COMMON/PLOTS/NDIM,XMIN,YMIN,DX,DY,NBINX,NBINY,TITLE
	REAL*4 RS,RV,RQ,RH
	CHARACTER*4 RRS,RRV,RRQ,RRH
	EQUIVALENCE(RS,RRS),(RV,RRV),(RQ,RRQ),(RH,RRH)
	DATA RRS/'   S'/
	DATA RRV/'   V'/
	DATA RRQ/'   Q'/
	DATA RRH/'   H'/
100	NDLTY=1000
	NPLTS=1
C	CALL ASSIGN(4,'TI:')
C	CALL ASSIGN(5,'TI:')
C	CALL ASSIGN(6,'LP:')
	IXTR=0
	DO 1982 N=1,128
1982	LINE(N)=0
	WRITE(0,8000)
8000	FORMAT(' Give Output Dataset Name>',\)
	read(0,2)ionm
	DO 8222 N=1,50
	NNN=51-N
	IF(ICHAR(IONM(NNN)).GT.32)GOTO 8223
	IONM(NNN)=0
8222	CONTINUE
8223	CONTINUE
C	IDL=NNN+1
	CALL WASSIG(6,IONM)
	ITTFG=0
	IF((IONM(1).EQ.'C'.OR.IONM(1).EQ.'c').AND.IONM(4).EQ.':')
     1  ITTFG=1
C ALSO MAKE TTY IMAGES IF 1ST 2 CHARS ARE TT
	IF(IONM(1).EQ.'T'.AND.IONM(2).EQ.'T')ITTFG=1
	IF(ITTFG.EQ.1)IXTR=1
C TOGGLE FOR PLOT ROUTINE...
	WRITE(0,1)
1	FORMAT(' Enter plot command>',\)
	READ(0,2)LINE
2	FORMAT(128A1)
C FIND END OF LINE ENTERED BY LOOKING FOR 1ST CHAR BIGGER THAN SPACE IN ASCII.
	DO 1980 N=1,128
	NN=129-N
	IF(ICHAR(LINE(NN)).GT.32)GOTO 1981
1980	CONTINUE
1981	LQ=NN+1
	LOGF1=0
	LOGF2=0
	NBFG1=0
	IF(LINE(1).EQ.'P')NBFG1=1
C NBFG1 MAKES YMIN=0. THUS IF CMD STARTS WITH PP PLOT IS POSITIVE
C DITTO NBFG2
	NBFG2=0
	IF(LINE(2).EQ.'P')NBFG2=1
	IF(LINE(1).EQ.'L')LOGF1=1
	IF(LINE(2).EQ.'L')LOGF2=1
	LLA=LINE(1)
	LLB=LINE(2)
C 1ST 2 CHARS SAY LOG OR LOGLOG (IF 2DIM GRAPH)
C LOGF1 WILL TAKE LOG OF VEC1 AND LOGF2 WILL TAKE LOG OF VEC2
C IF SET.
C NOTE THAT THIS ALSO TAKES ABS OF NUMBER.
	LQ=LQ+1
	LQ=MIN0(128,LQ)
	LINE(127)=0
	LINE(128)=0
	LINE(LQ)=0
c
c process switches.
c switches are after trailing + sign
c
c +hnnn = set height
c +wnnn = set width
	nhov=0
	nwov=0
	KK=INDEX(LINE,'+')
	IF(KK.GT.50)GOTO 6000
C SKIP THIS AREA IF NO SWITCHES ARE FLAGGED
	LINE(KK)=CHAR(0)
C SKIP SWITCHES IN LATER PROCESSING.
c since we look for a number, first try to decode the number as a
c 3 digit one...
	kkk=kk+2
	lend1=kkk+30
	call gn(kkk,lend1,num1,line)
c num1 can be h or w depending on line(kk+1)
	if(line(kk+1).eq.'h'.or.line(kk+1).eq.'H')nhov=num1
	if(line(kk+1).eq.'w'.or.line(kk+1).eq.'W')nwov=num1
C GN RETURNS ITS LAST CHAR AFTER THE # IN ITS 1ST ARG.
	IKK=INDEX(LINE(KK+1),'+')
	IF(IKK.GT.30)GOTO 6000
	KKK=IKK+KK
c 2nd + sign flags 2nd switch...
	kk=kkk+2
	lend1=kk+30
	call gn(kk,lend1,num1,line)
	if(line(kkk+1).eq.'h'.or.line(kkk+1).eq.'H')nhov=num1
	if(line(kkk+1).eq.'w'.or.line(kkk+1).eq.'W')nwov=num1
c that should do it...
6000	CONTINUE
	LS1=INDEX(LINE,32)
C CALL OUR PORTACALC INDEX FCN
	KK=LS1+1
	LS2=INDEX(LINE(KK),32)
	IF(LS1.GT.40.OR.LS2.GT.40)WRITE(0,25)LS1,LS2,LQ
25	FORMAT(' Spaces not seen. Find spaces at ',3I6,
     1  /,' Usage: PCG file V1:V2 C V3:V4 +HNNN+WNNN')
	IF (LS1.GT.40.OR.LS2.GT.40)GOTO 100
	LINE(LS2+LS1)=0
	CALL RASSIG(1,LINE(LS1+1))
C SET UP FILE 1 TO READ SAVED FILE FROM SHEET
	LINE(LS2+LS1)=32
	LX=LS1+LS2+1
C SCAN THE REST STARTING AT LX
C GRAB OFF OUR ARGUMENTS FIRST, THEN GET ON WITH THE PLOTS.
	CALL PLOT(0.,0.,-1,0)
C HOWEVER INITIALIZE PLOT ARRAY EARLY ON.
	K1=LX
	K2=110
	CALL GVSCAN(LINE,K1,K2,LSTCHR,ID1,ID2,IVLD)
	IF (IVLD.NE.0)GOTO 150
	WRITE(0,3)
3	FORMAT(' First variable invalid. Try again.')
	GOTO 100
150	CONTINUE
	IF(LINE(LSTCHR).EQ.':')GOTO 160
	WRITE(0,4)
4	FORMAT(' Colon missing in first range.')
	GOTO 100
160	CONTINUE
	K1=LSTCHR+1
	K2=110
	CALL GVSCAN(LINE,K1,K2,LSTCR,ID1B,ID2B,IVLD)
	IF (IVLD.NE.0)GOTO 164
	WRITE(0,5)
5	FORMAT(' 2nd variable in 1st range invalid.')
	GOTO 100
164	CONTINUE
	IF(ID1.NE.ID1B.AND.ID2.NE.ID2B)GOTO 166
	GOTO 167
166	WRITE(0,6)
6	FORMAT(' Variable pair not in a row or column together')
	GOTO 100
167	CONTINUE
	KCR=1
C : CHECK FOR '' OLD VERSION..........
	IF(LINE(LSTCR).EQ.'[')GOTO 170
	LSTCR=LSTCR+1
	IF(LINE(LSTCR).EQ.'[')GOTO 170
169	WRITE(0,7)KCR
7	FORMAT(' Invalid format of [c] character ',I5)
	GOTO 100
170	LSTCR=LSTCR+1
	KCR=2
	IF(LINE(LSTCR).EQ.']')GOTO 169
	KLET=LINE(LSTCR)
	LSTCR=LSTCR+1
C SCAN OVER NEXT ']' NOW
	KCR=3
	IF(LINE(LSTCR).NE.']')GOTO 169
	LSTCR=LSTCR+1
C IF WE PICK UP A VALID VARIABLE HERE, ALL'S WELL. OTHERWISE WE HAVE
C A HISTOGRAM AND WE'RE DONE (FOR THIS VERSION ANYHOW)
	K1=LSTCR
	K2=110
	NDIM=1
	CALL GVSCAN(LINE,K1,K2,LSTT,ID1C,ID2C,IVLD)
	IF(IVLD.EQ.0)GOTO 200
C IF HERE, THERE HAS TO BE 1 MORE VARIABLE DECODED AND TESTED.
	IF(LINE(LSTT).EQ.':')GOTO 175
	WRITE(0,8)
8	FORMAT(' Invalid second variable range.')
	GOTO 100
175	CONTINUE
	K1=LSTT+1
	K2=110
	CALL GVSCAN(LINE,K1,K2,LSTCC,ID1D,ID2D,IVLD)
	IF(IVLD.NE.0)GOTO 180
	WRITE(0,9)
9	FORMAT(' Invalid 2nd variable of 2nd range')
	GOTO 100
180	CONTINUE
C NOW ALL DECODED.
	NDIM=2
C NOW WE HAVE SET UP THE DIMENSION OF OUR PLOT.
200	CONTINUE
C NOW IT'S POSSIBLE TO READ IN THE FILE ONCE TO NORMALIZE IT, THEN
C REWIND AND READ AGAIN TO PLOT IT.
	XMIN=99.E10
	YMIN=99.E10
	IF(NBFG1.NE.0)YMIN=0.
	IF(NBFG2.NE.0)XMIN=0.
C SET TERRIBLY LARGE X,Y MINS UNLESS POSITIVE PLOT, THEN START AT 0.
C (WE'LL FIX THEM UP!)
	XMAX=-99.E10
	YMAX=-99.E10
C SET UP MAXIMA ALSO IN BOGUS WAY. THIS ENSURES WHATEVER WE GET
C WILL BE BETTER THAN OUR "FIRST GUESS".
C  INSERT TITLE AS OUR COMMAND LINE, FOR INTERNAL DOCUMENTATION.
C (LX THRU END)
	DO 11 N=1,78
11	LTITL(N)=CHAR(32)
	LX=LS1
C INCLUDE FILENAME TOO.
	DO 12 N=1,50
	LTITL(N)=LINE(LX)
	IF(LX.GT.76)GOTO 13
12	LX=LX+1
C FLAG LOG SCALE FLAGS IN TITLE
	IF (LLA.EQ.'L')TITLE(18)='LOGX'
	IF (LLB.EQ.'L')TITLE(19)='LOGY'
13	READ(1,10)LINE
10	FORMAT(128A1)
	IF(NDIM.EQ.2)GOTO 17
	XMIN=0.
	XMIN2=0.
	ICNT=0
17	CONTINUE
C IGNORE TITLE, JUST READ IT IN, THEN FORGET IT.
	IV1=1
	IV2=1
220	CONTINUE
C NOTE WE READ IN THE NUMBER IN NUMERIC FORMAT. EASIER THAT WAY.
	irrw=0
	iccl=0
	READ(1,14,END=250,ERR=224)LET1,IRRW,ICCL,XYVAL
224	continue
14	FORMAT(A1,I5,1X,I5,1X,E50.35)
	READ(1,15,END=250,ERR=225)LFVLD,(LINE(IV),IV=120,128),KKTYP
225	continue
15	FORMAT(I3,1X,9A1,1X,I5)
C READS IN AN ENTRY OF SAVED SHEET. TEST IF IN OUR RANGE.
	IF(IRRW.GE.ID1.AND.IRRW.LE.ID1B.AND.ICCL.GE.ID2.AND.ICCL
     1  .LE.ID2B)GOTO 221
	IF(NDIM.NE.2)GOTO 223
	IF(IRRW.GE.ID1C.AND.IRRW.LE.ID1D.AND.ICCL.GE.ID2C
     1  .AND.ICCL.LE.ID2D)GOTO 222
	GOTO 223
221	CONTINUE
C NUMBER IS IN FIRST RANGE TO PLOT. FIGURE IT OUT.
	IF(LOGF1.NE.0.AND.XYVAL.NE.0)XYVAL=ALOG10(ABS(XYVAL))
	VEC1(IV1)=XYVAL
	IV1=IV1+1
	IF(NDIM.EQ.1)ICNT=ICNT+1
	IF(NDIM.EQ.1)XMAX=ICNT
	IF(NDIM.EQ.1)GOTO 18
	IF(XYVAL.LT.XMIN)XMIN=XYVAL
	IF(XYVAL.GT.XMAX)XMAX=XYVAL
	GOTO 223
18	CONTINUE
	IF(XYVAL.LT.YMIN)YMIN=XYVAL
	IF(XYVAL.GT.YMAX)YMAX=XYVAL
	VEC2(IV2)=FLOAT(ICNT)
	IV2=IV2+1
	GOTO 223
222	CONTINUE
	IF(NDIM.EQ.1)GOTO 223
	IF(LOGF2.NE.0.AND.XYVAL.NE.0)XYVAL=ALOG10(ABS(XYVAL))
C NUMBER IS IN SECOND RANGE SELECTED.
C KNOW IT'S A Y COORD HERE.
	VEC2(IV2)=XYVAL
	IV2=IV2+1
	IF(XYVAL.LT.YMIN)YMIN=XYVAL
	IF(XYVAL.GT.YMAX)YMAX=XYVAL
223	CONTINUE
	GOTO 220
250	CONTINUE
C NOW MINIMA,MAXIMA ALL SET UP.
	IF(XMIN.GE.XMAX.OR.YMIN.GE.YMAX)CALL ENDXQQ
C EXIT IF NOTHING IS THERE TO GRAPH.
	XRANGE=XMAX-XMIN
	YRANGE=YMAX-YMIN
C
	IF(XRANGE.LE.0)XRANGE=60.
	IF(YRANGE.LE.0)YRANGE=20.
C	XNUM=100.
	AMXRG=100.
	AMYRG=50.
	IF(ITTFG.EQ.1)AMXRG=60.
	IF(ITTFG.EQ.1)AMYRG=20.
	XNUM=AMXRG
	IF(NDIM.EQ.1.AND.(XRANGE.LT.100.))XNUM=XRANGE
	YNUM=AMYRG
	IF(NDIM.EQ.2)XNUM=AMYRG
C HANDLE SWITCHES THAT OVERRIDE HEIGHT AND WIDTH TO USE.
	IF(NHOV.NE.0)YNUM=NHOV
	IF(NWOV.NE.0)XNUM=NWOV
	DX=XRANGE/XNUM
	DY=YRANGE/YNUM
	IF(.NOT.(NDIM.EQ.1.AND.DX.LT.1))GOTO 19
	IF(NWOV.EQ.0)DX=1.
19	NBINX=XNUM
	NBINY=YNUM
	CALL PLOT(RV,0.,0,1)
C INITIALIZE PLOT
C NDIM, MINIMA, MAXIMA ALL SET UP NOW.
C
C WE SAVED VALUES IN VEC1,VEC2 AND PLOT THAT WAY.
C  ALSO NOTE BOTH ALWAYS EXIST.
C
	LENGTH=MIN0(IV1,IV2)
C SAME IF NDIM=1
	DO 20 N=1,LENGTH
	IF(NDIM.EQ.1)CALL PLOT(VEC2(N),VEC1(N),1,1)
	IF(NDIM.NE.1)CALL PLOT(VEC1(N),VEC2(N),1,1)
20	CONTINUE
C PLOT IT OUT NOW
C CHOOSE OPTION FOR FORMAT (SCALE, VARY HEIGHT, SHADE)
	X=RS
	IF(KLET.EQ.'V')X=RV
	IF(KLET.EQ.'Q')X=RQ
	IF(KLET.EQ.'H')X=RH
	CALL PLOT(X,0,2,1)
	CALL ENDXQQ
	END


