)	SUBROUTINE BINP2 (K,IX,IY,IZ,NA,ISUM,NP):C_Title  BINP2 bin, histogram and display for  I*2 arrays.C_Arguments8C All arguments below are dependent on the calling mode.7	INTEGER K	![I] Varies with mode see description below.6	INTEGER*2 IX(*)	![I] Abscissa of values to be binned.6	INTEGER*2 IY(*) ![I] Ordinate of values to be binned.2	INTEGER*2 IZ(*) ![I] Values to be summed in bins.3	INTEGER*2 NA(*) ![I] Number of points in each bin.3	INTEGER*4 ISUM(*) ![B] Sum of values in each bin. <	INTEGER NP	  ![O] Number of points in the  IX,IY,IZ arrays.GC  Argument meanings depend on the value of mode (first argument). ****NC  They are most representative for mode 3, where they are in compilable form,+C  but are described here for all 6 modes. JC  Where the variable name is used in this explanation, the same variable 'C  must be used in the calling program.CC Mode 1, initialization.%CC K equal to 1, Initialization mode.)CC IXLIM(5)   X range and number of bins.0CC	     1 = Minimum value of  X, in print units.0CC	     2 = Maximum value of  X, in print units.1CC	     3 = NX is equal to the number of  X bins.%CC	        (not counting the border)./CC	     4 = Minimum value of  X, in disk units./CC	     5 = Maximum value of  X, in disk units.=CC	IYLIM(5) Input, same as the input items for X, but for  Y.9CC	IOP	 Input, logical unit for print file (when mode=5).BCC	NA(*)    Output, count array. will be zeroed. Must be at least :CC       NX*NY   no border, or (nx+2)*(ny+2) for a border.FCC	ISUM(*)	Output, sum or average the array, this array may be zeroed.&CC		same size requirements as for  NA.DCC 	        ISUM will not be zeroed if it is the same array as  NA. @CC		thus, if only the count will be constructed, and calls with NCC               mode equal to 3 or 4 are not planned then this can be called 5CC 		with the  NA array which is also used for  ISUM.7CC		and there is no requirement for  INTERGER*4 arrays.DCC	LP2	Input border flag, 0 indicates no border, otherwise a border.CCC Mode 2, add to the count.,CC	K	Input mode indicating add to the count./CC	IX(*)	Input abscissa of values to be binned.1CC	IY(*)   Input ordinate of values to be binned.CC	IDUMMY 	Not used.,CC	NA(*) 	Both number of points in each bin.CC	IDUMMY  Not used.2CC	NP	Input number of points in the  IX,IY arrays.C%CC Mode 3, add to the count and sums.#CC	K	Input, mode to count and sums./CC	IX(*)	Input abscissa of values to be binned.0CC	IY(*) 	Input ordinate of values to be binned.,CC	IZ(*) 	Input values to be summed in bins.-CC	NA(*) 	Input number of points in each bin.3CC	ISUM(*) Input,Output sum of values in each bin. 5CC	NP	Input number of points in the  IX,IY,IZ arrays.CCC Mode 4, average the sums."CC	K	Input mode. Average the sums.CC	IDUMMY  Not used.HCC	IEXCL	Input value to be used as the average if the bin count equal 0.'CC	IZ(*)   Output average value in bin..CC	NA(*)   Input number of points in each bin.,CC	ISUM(*) Input sum of values in each bin. <CC	IFAC	Input divisor to be used in scaling averaged values.CC		0 = no scalingCC Mode 5, print the results.'CC	K	Input, 5 equals print the results."CC	BUF(120) Output scratch buffer.=CC	IREP(10) Input, Output control representation; and report.:CC		(1) Input 0 results in write in normal integer format.*CC			 -1 blank all occurences of  IREP(2)..CC			 greater than 0 make a crude shade print.-CC                        (see  L in  LDENS).;CC		(2) Input, exclude value, or the minimum value for the !CC                   shade print.CC		(3) Not used@CC		(4) Input constant to be added to the  Z values for display.=CC		(5) Input modulo for  Z, or (if irep(1) greater than 00) (CC                   shadeing step size.0CC		(9,10) Output total number of points binned.+CC	IZ(*) 	Input average values in each bin.1CC	MTITLE(20) Input 40-byte title of the dataset.ICC	IFORM   Input bin integer format length (1 to 6), 0 equals no display.?CC	ISP	Input number of printer lines to space between each row.ACC		If ISP is 0=don't space,  1 to 9 space n lines between print,<CC		and add 10 to print just the  Y value for a line of all HCC 	        empty bins; Add 20 to print nothing (not even the  Y value) CC		for empty lines.CC Mode 6, remove the border.7CC	K	Input mode, 6 operation code to remove the border.ECC	IBOR(*) Input bordered array (NX+2)*(NY+2), could be either count CC               or average.8CC	IUNB(*) Output array after removal of border (NX*NY).CC	IDUMMY  Not used.CC	IDUMMY  Not used.CC	IDUMMY  Not used.CC	IDUMMY  Not used.
C_Description<C This is a six-part routine.  ARGUMENTS  VARY  WITH  MODE !JC This routine is meant to provide all aspects of constructing and simple KC displaying of 2 or 3 dimensional data.  This routine also allows scaling &C of the printout for all 3 variables.?C The standard usage is to call this routine, for each dataset:FC Call once with  K equal to 1 to initalize. Call any number of times )C with  K=2 or  K=3 to add data points.  ?C Call BINP2 with mode 4 (K equal to 4) to average in the bins.IC Call one or more times with  K=5 to display the data (several options),EC and call once with  K=6 to remove a border before saving the array.LC The interior to bins is greater than the minimum and less than or equal toNC the maximum.  i.e., Upper bin boundaries are inclusive and lower boundaries C are exclusive.JC Y-axis labeling: Value refers to lower boundary of that line. therefore,6C value is the maximum included in the next lower bin.LC X_axis labelling (every 12 character locations). Value is for the boundary)C immediately to the right of the + sign.>C Warning: This routine will round down non-integer bin width!(C_Call  A2AN  FILL  LDENS  LIMIT2  ZSUP22C_Hist	78feb20  HUGH_H_KIEFFER  U.S.G.S._FLAGSTAFFC	81apr08  HHKGC	83jan06  R.MEHLMAN  FOR  XBIN with minor changes and much respelling C	         and KKI output.IC	87oct14  HHK revise documentation to  NIMS/VIMS, and greatly elaborate..C	         also, start toward structured code.@C	87nov13  HHK revise  FORMAT construction, and number of  TICs.C_End	INTEGER*4 KK,KKI 	INTEGER*2 IEXCL,IZMOD,IZCON, NFF	INTEGER*2 IXX(11),IW(122)      !Buffers for x-values and output line.,	INTEGER*2 IDUM2(3) 	       !Scratch buffer.?	BYTE LW(132)		       !Character buffer for  ascii output line.C	LOGICAL LP2	!True = Add a border around array for external points.7	LOGICAL LSUM	!True = Process values, as well as count.$C constructed format for printout.  3C          1 2 3 4 5 6 7 8 9 101112     1314151617 3C  IFMT = "( 1X,I5,I2,'+', _NF_ I(IFORM),'+', I2) "*C		    2X           A 1                   <	INTEGER*2 IFMT(16)/2H(  ,2H1X ,2H,I ,2H5, ,2HI2 ,2H,' ,2H+';     &,2H,  ,2H   ,2H  ,2HI ,2H 6 ,2H,' ,2H+' ,2H,I ,2H2) /2C insert for field length. IFMW = 1  2  3  4  5  66	INTEGER*2 IFMW(6)/2H 1 ,2H 2 ,2H 3 ,2H 4 ,2H 5 ,2H 6/"	CHARACTER*12 TIC /'+...........'/	CHARACTER*1 PLUS /'+'/!C  Branch on mode dependent on K.#	GO TO (100,200,300,400,500,600),K	9C _1_setup (1,xlim(5),ylim(5),printer_lun,count, sum,lp2)8CC   BINP2 (K,     IX,     IY,	       IZ,   NA,ISUM,NP );C limits names:  I=disk,  M=print units  D=delta,  O=originFC for no border,loc=((x-min)+(d-1))/d. so orig=min+1-d [-d for border]0C     IDX = delta x unit for bins in disk units.0C     IDY = delta y unit for bins in disk units.&C     IOX = origin of x in disk units.&C     IOY = origin of y in disk units.&C     MDX = Delta of x in print units.&C     MDY = Delta of y in print units.C     IOP = printer units.C     NX = Number of X bins.100	IOP=IZ(1)	        	NX=IX(3)			NXP2=NX                  	IDX=(IX(5)-IX(4))/NX	IF (IDX.EQ.0) IDX = 1	IOX=IX(4)+1-IDX	MDX=(IX(2)-IX(1))/NX
	MOX=IX(1)2C-------Do the same for  Y, NY = number of y bins.		NY=IY(3)	NYP2=NY	IDY= (IY(5)-IY(4))/NY     	IF (IDY.EQ.0) IDY = 1	IOY=IY(4)+1-IDY	MDY=(IY(2)-IY(1))/NY
	MOY=IY(1)0C  JT1 = index of first y bin not in the border.0C  JT2 = index of last  y bin not in the border.?C  LP2 = Flag for border, set true if input value was not zero.	JT1=1				JT2=NY			LP2 = (NP.NE.0)	9C    NXP2 = NX bins plus 2, if a border around main grid.9C    NYP2 = NY bins plus 2, if a border around main grid..C    JT2 = Number of intervals in Y direction.	IF (LP2) THEN		  NXP2=NX+2	  NYP2=NY+2	  IOX=IOX-IDX          IOY=IOY-IDY	  MOY=MOY-MDY          JT1=2
	  JT2=NYP2-1	ENDIF/	NXY=NXP2*NYP2           !Total number of bins.@	ISUM(1)=-1	        !Set a test value in the first word of ISUM.,	CALL FILL (NA,NXY,0)	!Zero the count array.7	ITEST=ISUM(1)		!unless first word of ISUM is now zero.C    Zero sum array.+	IF (ITEST.NE.0) CALL FILL(ISUM,NXY+NXY,0)		KK=0	KC=0              .C    KKI = Number of points within the border.	KKI=0	                ?C     Set up limit2.  "FLAGONLY" is true if there is no border.,	CALL LIMIT2 (.TRUE.,1,NXP2,1,NYP2,.NOT.LP2)	GO TO 9HC  KC=#calls  KK=#points  KKI=#points in bounds  (=KK for bordered plot)/C _2_add to count  BINP2 (2,x,y,-,count, - ,#p)200	LSUM= .FALSE.
	GO TO 310,C _3_ &   to sum  (3, x, y, z,count, sum,#p)&C	.  BINP2  (K,IX,IY,IZ,   na,isum,np)300	LSUM=.TRUE.310	KC=KC+1	IF(NP.LE.0) GO TO 9AC     LC = x and y coordinates relative to origin in print units.DC     determine if point is in the bounds of our grid, if so convert!C     to grid location and count.
	DO I=1,NP          LCTX=(IX(I)-IOX)/IDX	  LCTY=(IY(I)-IOY)/IDY6          CALL LIMIT2 (.FALSE.,LCTX,LCTY,ITX,JTY,LIM2)          IF (LIM2) THEN	    KKI=KKI+1H	    J=(LCTY-1) * NXP2 + LCTX            !Convert x_y location to index.1	    NA(J)=NA(J)+1	                !Add to count.5	    IF (LSUM) ISUM(J) = ISUM(J) + IZ(I)	!Add to sum.		   ENDIF         ENDDO		KK=KK+NP	IF (LP2) KKI=KK	GO TO 90C _4_average(4, -,iexcl,average,count, sum,ifac),C BINP2 (K,IX,   IY,     IZ,   NA,ISUM,  NP))C sets average to iexcl if count is zero.C Set IFAC as z scaling factor.400	IFAC=1			IF (NP.NE.0)IFAC=NP	IEXCL=IY(1)	DO J=1,NXY          IF(NA(J).NE.0) THEN 	     IZ(J)=ISUM(J)/(NA(J)*IFAC)          ELSE	     IZ(J)=IEXCL          ENDIF
        ENDDO	GO TO 9.C _5_plot(5,buf,iy1-5,iz,mtitle(20),iform,isp).C  BINP2 (k, ix,   iy,iz,        na, isum, np)$C LOUT = 10's digit: 0 = full print.FC	             1 = Y value only for empty bins, 2 = bypass empty line.5C LISP = 1's digit = lines to space between printing.500	LOUT = NP / 10		LISP = NP - 10 * LOUT	CC	IF (KKI.GT.0) WRITE(IOP,501)501	FORMAT('1');C-------Write upper header, return number of points binned.)	WRITE (IOP,301) KC,KK,KKI,(NA(I),I=1,30)A301	FORMAT(' BINP2:',I6,'=CALLS ',I8,'=POINTS ',I8,'=IN',4X,30A2)	CALL A2AN(KKI,IY(9),2)	IF (KKI.LT.1) GO TO 9IC-------Setup for printout. transfer arguments in  IY to local variables..C       Format of 0 means bypass all printing.	CALL A2AN (ISUM,IFORM,1)	IF (IFORM.LE.0) GO TO 9	     	IF (IFORM.GT.6) IFORM=6
	LBLANK=IY(1)	IEXCL=IY(2)	IZCON=IY(4)	IZMOD=IY(5)		IS=IZMOD!C  No modulo for density shading.	IF (LBLANK.GT.0) IZMOD=0	
	NXPP=NXP2	MAXP=120/IFORM +2	IF (NXPP.GT.MAXP) NXPP=MAXPLC  Setup format. only the interior bins may be converted into non-I formats.	IF (LBLANK.EQ.0) THEN)		IFMT(11) = ' I'		!Set output word type.&		IFMT(12) = IFMW(IFORM)	!field width.'		NF = NX			!number of interior fields.		    ELSE)		IFMT(11) = ' A'		!set output word type. 		IFMT(12) = '1 '		!field width.3		NF = IFORM*NX		!number of interior output fields.
	    ENDIF;	  CALL ASCII (NF,LW)		 !convert number of fields to ascii.>	  CALL MVBYTE (LW,3,4,IFMT(9),1) !move into format statement.&	IFMT(5) = '2X'			 !left exterior bin..	IF (LP2) IFMT(5) = 'I2'		 !left exterior bin.,C determine how many  X values to be printed	NTIC = (NX*IFORM+11)/12 + 1 	IF (NTIC.GT.11) NTIC = 11!C---    Convert to phsical units.	DO I=1,NTIC+	  IXX(I) = MOX + ((I-1) * 12 / IFORM) *MDX-	ENDDO                                       "	WRITE (IOP,303) (IXX(I),I=1,NTIC)303	FORMAT (2X,I8,10I12)304	FORMAT (8X,12A)!C------Print out the stored plot.OC  BIG loop on lines.__________________________________________________________5C  IXO = location in array of first bin in this line.+C  IXL = Location of last bin in this line.!C  IV = Value of Y for this line.	DO J=1,NYP2			  IXO = (NYP2-J) * NXP2 + 1		  IXL = IXO + NXP2 - 1		  IV = MOY + (NYP2-J) * MDY 	  IF (J .EQ. NYP2) IV=0C  Write upper axis.6	  IF (J.EQ.JT1) WRITE(IOP,304) (TIC,I=1,NTIC-1),PLUS	
          M=0+C  Will be right-most bin that is not empty
	  NK=0			%C  Transfer from bins to line buffer.AC  Find last filled bin IZ(I), add a constant to Z placing in IW.	 DO I=IXO,IXL           M=M+1	   IF (IZ(I).NE.IEXCL) NK=M		   IW(M)=IZ(I)+IZCON	 ENDDO%C  Truncate left digits if requested.	IF (IZMOD.NE.0) THEN		    DO I=1,NXP2	      IW(I)=MOD(IW(I),IZMOD)
	    ENDDO	 ENDIFC  Insert blank lines.	IF (LISP.GT.0) THEN		  DO I=1,LISP	    WRITE (IOP,303)	  ENDDO	ENDIF'C  If no short flag, set to full-right 	IF (LOUT.EQ.0) NK=NXP2		IF (NK.EQ.0) THEN	/C   All bins in this line are empty (excluded).%	   IF (LOUT.LT.2) WRITE(IOP,IFMT)IV	&C  LBLANK = flag for denstity shading.*C  Convert to  ASCII for zero-suppression.+C  NOUT = number of ascii locations output.EC  NKK = number of interior  bins to be converted to density shading.	  ELSEIF (LBLANK.NE.0) THEN 	    NKK = NK + 1 - JT1			    IF (NKK.GT.NX) NKK = NX		    IF (NKK.LT.1) NKK = 1	    NOUT = NKK * IFORM		+C  Zero (i.e., exclude value) -suppression.	    IF (LBLANK.LT.0) THEN :	      CALL ZSUP2 (IW(JT1),NKK,IEXCL,IFORM, LW,IDUM,IDUM2)
	    ELSE	'C  Density shading.  NA is MTITLE here.;	     CALL LDENS (IW(JT1),NKK,IFORM,LBLANK,NA,IEXCL,IS, LW)
	    ENDIFGC  LP2 true if possiblility of borders, if false there are no borders,  C  write all data that is valid."C  If NKK = NX then interior full.1C  else valid data does not extend to full right.	   IF (LP2) THEN 		      IF (NKK.EQ.NX) THEN	6		 WRITE (IOP,IFMT) IV,IW(1),(LW(I),I=1,NOUT),IW(NXP2)
	      ELSE		-		 WRITE (IOP,IFMT) IV,IW(1),(LW(I),I=1,NOUT)	      ENDIF
	   ELSE		+	      WRITE (IOP,IFMT) IV,(LW(I),I=1,NOUT)		   ENDIF	  ELSEC  Output in integer format.'	    WRITE (IOP,IFMT) IV,(IW(I),I=1,NK)	  ENDIFC  Write lower axis.4	IF (J.EQ.JT2) WRITE(IOP,304) (TIC,I=1,NTIC-1),PLUS		ENDDO@C  End of line loop.____________________________________________C  Write lower header."	WRITE (IOP,303) (IXX(I),I=1,NTIC)	GO TO 98C _6_remove border (6,bordered,unbordered,- , -,   -, -)-C	BINP2 (k,ix,      iy,        iz,na,isum,np)FC LP2 is flag for the border, if false don't bother because the borderC must be there to remove!600	IF (.NOT.LP2) GOTO 9	3C Remove the border, keeping only the data in bins.
	DO I=1,NY	   I1 = 2 + I * NXP2	   I2 = 1 + (I-1) * NX 	   CALL A2AN (IX(I1),IY(I2),NX)	ENDDO9	RETURN	END