	PROGRAM SPREP4
C
C	Version 5.0, October 8, 1997, by M. Knappertsbusch
C
C	Preparation for shape analysis.
C	!!! Processing in batch mode. !!!
C
C	Modifications:
C	19.12.2000: Addition of error reports for the cases if 1.) files are missing,
C	            2.) if the number of points of the input files are less than
C	            the number of new points (=M), and 3.) if the number of interpolated
C	            points written is not identical to M.
C
C	2.1.2000:   Reorganization and formatting of output file MEASUREMENTS (omit output of
C	            XMIN,XMAX,YMIN,YMAX, XCTR,YCTR,DIST,DEGREES,DEGREES and only output of
C	            FILENAME, Age, N, XMAX-XMIN, YMAX-YMIN, and AREA. Calculation of DIST is
C	            outcommented.
C
C	6.3.2001:   Minor modification around error message if the # of outlines is less than
C		    the new # of points (250). [See around statement labels 10,12 and 15].
C
C	This version interpolates points in anti-clockwise direction.
C	It can also handle X and Y coordinates in pixels or µm.
C	First calculate minima and maxima of X and Y coordinates
C	of the input outline, and calculate coordinates of centroid.
C	The X and Y coordinates are then transformed into polar coordinates
C	arount the centroid.
C
C	Input files are:
C		- File with name TRACED_FILES containing a list of all
C		  names of files, that contain the traced X,Y coordinates
C		  (in cartesian coordinates).
C		- All files, that contain the traced X,Y coordinates
C		  of the outlines. Remember: The filename convention is
C		  that the sample files consist of 15 characters, and that
C		  suffixes are separated by a dash.
C
C	Output files are:
C		- Files with the suffix _POL, which contain the polar coordinates
C		  of the outline. All angular values (THETA) are in Radians.
C		- Files with the suffix _INT contain a reduced set of X,Y outline points
C		  which are interpolated at equiangular distances. Values are in
C		  cartesian coordinates.
C		- A file called "Errors_during_Sprep" recording errors during processing.
C		- The file with the name "MEASUREMENTS" contains individual size and
C		  shape parameters for each file:
C
C		  FILENAME,AGE,N,XMAX-XMIN,YMAX-YMIN,Area
C
C	Calculations:
C     Outcommented. - Maximal diameter of the closed plane curve (DIST) is calculated.
C     - The area of the light object (AREA) is calculated.
C
C	Note, that the input file must be terminated by a carriage
C	return (otherwise the program disregards the last record).
C
C
	INTEGER N,M,NINT
	LOGICAL FILEX
	DOUBLE PRECISION XMIN,XMAX,YMIN,YMAX,XCTR,YCTR
	DOUBLE PRECISION X,Y,U,V,S,PHI
	DOUBLE PRECISION DELTA,THETN,RN,THETN1,RN1,R,THETA
C	DOUBLE PRECISION DIST,OMEGA1,OMEGA2,ANGLE1,ANGLE2	! Outcommented (max diameter is not determined)
	DOUBLE PRECISION AREA,ARC,R0,THET0
	DOUBLE PRECISION ALPHAN,RHON,EPS
	DOUBLE PRECISION PI
	REAL AGE
	CHARACTER*1  ANSWER
	CHARACTER*11 HEADER
	CHARACTER*20 FILE_LIST
	CHARACTER*17 INPUT
	CHARACTER*19 OUTPUT1,OUTPUT2
C
	PARAMETER(PI=3.14159265)
C
C
	WRITE(9,*) 'Preparation for outline analysis:'
	WRITE(9,*) '. . .Enter the age of the sample (Ma). . .'
	WRITE(9,*) '       (maximum 5 characters long)'
	READ(9,*) AGE
	WRITE(9,*) '. . .Enter new number of points (250). . .'
	READ(9,*) M
C
C	Open file with name "MEASUREMENTS" and write header to it:
C
	OPEN(14,FILE='MEASUREMENTS',STATUS='NEW')
	WRITE(14,2) AGE
 2	FORMAT(F5.2)
	WRITE(14,*) 'FILENAME, Age in Ma, # of points, XMAX-XMIN (µm),
     1 YMAX-YMIN (µm), Area (mm2)'
C
C	Recording errors during program run:
C
	OPEN(25,FILE='Errors_during_Sprep',STATUS='NEW')
C
C
C	Now, get filename from file list with name TRACED_FILES, open that file and
C	read the data:
C
	WRITE(9,*) '. . .Names of traced outlines needed:. . .'
	WRITE(9,*) 'enter list containing names of files (20 chars)'
	READ(9,4) FILE_LIST
 4	FORMAT(A20)
C
	OPEN(20,FILE=FILE_LIST,STATUS='OLD')
C	Determine sample name from INPUT, and write the sample name into the
C	header of file Errors_during_Sprep:
C
	READ(20,6) INPUT
	HEADER=INPUT(1:11)
	WRITE(25,3) HEADER
 3	FORMAT('Sample name: ',A11)
	REWIND(20)
C
C***	Main loop starts at statement number 5:
C
 5	READ(20,6,END=999) INPUT
	WRITE(9,7) '. . .Reading new file ',INPUT,'. . .'
 6	FORMAT(A17)
 7	FORMAT(A22,A17,A5)
C
C***	Check first, whether file with name INPUT exists or not. If it does not exist:
C	print a message and continue reading next file.
C
	INQUIRE(FILE=INPUT,EXIST=FILEX)
	IF (.NOT.FILEX) THEN
 	  WRITE(9,*) 'Error: file ',INPUT,' is missing'
	  WRITE(25,8) INPUT
	  WRITE(14,8) INPUT
 8	  FORMAT('File ',A17,' is missing')
	  WRITE(9,*) '. . .Program continues. . .'
 	  GOTO 5
	END IF
C
C
C	(Determination of actual filenames:)
C
	OUTPUT1=INPUT(1:15)//'_POL'
	OUTPUT2=INPUT(1:15)//'_INT'
C
C
C	Determination of size and shape parameters:
C
C*****	Determination of extrema:
C
	WRITE(9,*) '. . .Extrema are being calculated. . .'
C
	CALL EXTREMA(XMIN,YMIN,XMAX,YMAX,N,INPUT)
C	
C
C***	Record errors if the number of points (=N) in the Input file
C	is less than than the number of new points (=M, usually 250):
C	In case N<M the operator can choose whether this sample is to be included
C	in the analysis (normal continuation of the program) or whether this sample is
C	to be omitted and writing a message to file Errors_during_Sprep and continue
C	with next sample:
C
	IF (N.LT.M) THEN
	  WRITE(9,10) INPUT,N,M
	  WRITE(25,10) INPUT,N,M	! writing a message to file Errors_during_Sprep.	  
 10	  FORMAT('# Points in ',A17,' is ',I4,' i.e. < ',I4)
 	  WRITE(9,*) 'Include (=1) or skip and continue
     1 with next sample (=2) ?'
	  WRITE(9,*) 'Enter 1 or 2'
	  READ(9,12) ANSWER
 12	  FORMAT(A1)
 	  IF (ANSWER.EQ.'1') THEN
	    WRITE(9,*) '. . .Sample is included. . .'
	    GOTO 15
	  ELSE IF (ANSWER.EQ.'2') THEN
	    WRITE(14,10) INPUT,N,M	! writing a message to file MEASUREMENTS.
 	    WRITE(9,*) '. . .Sample skipped.
     1 Reading next sample. . .'
	    GOTO 5
	  END IF
	END IF
C
C
C*****	Dermination of X and Y coordinates of centroid:
C
 15	WRITE(9,*) '. . .Centroid is being calculated. . .'
C
	CALL CENTR(XCTR,YCTR,INPUT)
C
C*****	Transformation of X,Y coordinates into new cartesian 
C	coordinate system (U,V) with origin at XCTR,YCTR, and calculating
C	polar coordinates R AND THETA of U and V:
C
	OPEN(15,FILE=INPUT,STATUS='OLD')
	OPEN(16,FILE=OUTPUT1,STATUS='NEW')
C
 30	READ(15,*,END=500) X,Y
 	U=X-XCTR
	V=Y-YCTR
C
C	Now, calculation of polar coordinates of U and V and
C	write them to file "input name_POL":
C
	CALL POLAR(U,V,S,PHI)
	IF (PHI.LT.0.0) THEN
	  PHI=PHI+2.0*PI
	END IF
	WRITE(16,998) PHI,S
C***
	GOTO 30
 998	FORMAT(F12.8,',',F12.7)
C
 500	CLOSE(15)
	REWIND(16)
C
C
C*******Outcommented****************************************************************************
C	Calculation of the maximal diameter of the outline:
C
C	WRITE(9,*) '. . .Calculating maximum diameter. . .'
C
C	CALL DIAM(DIST,OMEGA1,OMEGA2)
C
C	Conversion into degrees:
C	ANGLE1=OMEGA1*180/PI
C	ANGLE2=OMEGA2*180/PI
C
C	REWIND(16)
C**********************************************************************************************
C
C
C	Now, calculation of the area of the light object: 
C	Coordinates are read from file with the polar coordinates. Variables used are
C	RN,RN1,THETN, and THETN1; these same variables are later superseeded
C	and re-used for the interpolation of points of the outline (directly
C	after the calculation of the area).
C
	WRITE(9,*) '. . .Calculate area of light object. . .'
C
C	First initialize AREA, and then read data from "input name_POL"
C	(THETN0,THETN, and THEN1 in Radians):
C
	AREA=0.0
	READ(16,*) THETN,RN
C
C	Store initial values separately:
C
	R0=RN
	THET0=THETN
C
C	Reading next values, calculating individual areas and summing up:
C
 40	READ(16,*,END=41) THETN1,RN1
	IF ((THETN.GT.1.5*PI).AND.(THETN1.LT.PI/2)) THEN
	  ARC=2*PI-THETN+THETN1
	ELSE
	  ARC=THETN1-THETN
	END IF
 	AREA=AREA+0.125*ARC*(RN1+RN)**2
	RN=RN1
	THETN=THETN1
	GOTO 40
C
C	Calculate last area and sum up:
C
 41	AREA=AREA+0.125*ARC*(RN1+RN)**2
C
	REWIND(16)
C
C	
C	Writing results to "MEASUREMENTS":
C
	WRITE(14,43) INPUT,AGE,N,XMAX-XMIN,YMAX-YMIN,AREA*1.0E-6
 43	FORMAT(A17,',',F5.2,',',I4,',',F6.1,',',F6.1,',',F7.4)
C
C
C
C	Now interpolation of R and THETA (in anti-clockwise direction):
C	All R,THETA,U and V are now superseeded by new values.
C	Given the succeeding polar coordinates RN,THETN and RN1,THETN1 
C	of a closed curve, this program interpolates along the curve, so
C	that the new points are at equiangular distances. The interpolated
C	points are written to file 'input name_INT'. All angular values are 
C	in Radians.
C
C	M is the new number of points along the curve.
C	DELTA is the angular spacing between new points.
C	THETA is the computed angular value for the next interpolated point.
C	U and V are the cartesian coordinates of the interpolated points
C	along the curve.
C
C	New values are written to file "input name_INT" (in cartesian 
C	coordinates U,V)
C
C
	WRITE(9,*) '. . .Interpolation at equiangular distances. . .'
C
	OPEN(17,FILE=OUTPUT2,STATUS='NEW')
C
C	Calculation of equiangular distance (=DELTA) from number of new points (=M):
C
	DELTA=2*PI/M
C
C	Read the starting point:
C
	READ(16,*) THETN,RN
C
C	Memorize first point:
C
	ALPHAN=THETN
	RHON=RN
C
	CALL IPOLAR(RN,THETN,U,V)
	WRITE(17,*) U,',',V
C
C	Count number of points (=NINT) written to OUTPUT2:
C	(This is for control purposes only, i.e. to check whether all
C	points from NINT=1 to M were written to OUTPUT)
C
	NINT=1
C
C	Initialize algoritm:
C
	THETA=THETN
 50	THETA=THETA+DELTA
C
C	Identity of 0 degrees and 360 degrees:
C
 	IF (THETA.GT.2*PI) THEN
	  THETA=THETA-2*PI
	END IF
C
C	Search in which new interval THETA is located and interpolate:
C	First, read next point:
C
 60	READ(16,*,END=600) THETN1,RN1
C
C	Consider singularity at transition 360 degrees --> 0 degrees:
C
	IF (((0.LT.THETN1).AND.(THETN1.LT.PI/2)).
     *    AND.(((3*PI/2).LT.THETN).
     *    AND.(THETN.LT.2*PI))) THEN
 70       CALL ROTATE(RN,THETN,RN1,THETN1,THETA,R)
	  CALL IPOLAR(R,THETA,U,V)
	  WRITE(17,*) U,',',V
	  NINT=NINT+1
	  EPS=THETA+DELTA
	  IF (EPS.GT.2*PI) THEN
	    EPS=EPS-2*PI
	  END IF
	  IF (((0.LE.EPS).AND.(EPS.LT.THETN1)).
     *      OR.((THETN.LE.EPS).
     *      AND.(EPS.LT.2*PI))) THEN
            THETA=THETA+DELTA
	    IF (THETA.GT.2*PI) THEN
	      THETA=THETA-2*PI
	    END IF
	    GOTO 70
	  END IF
	  RN=RN1
	  THETN=THETN1
	  GOTO 50
	END IF
C
C	All other points:
C
	IF ((THETN1.GE.THETA).AND.(THETA.GT.THETN)) THEN
	  CALL INTERP(RN,THETN,RN1,THETN1,THETA,R)
	  CALL IPOLAR(R,THETA,U,V)
	  WRITE(17,*) U,',',V
	  NINT=NINT+1
 80	  IF (THETN1.GE.ABS(THETA+DELTA).
     *    AND.ABS(THETA+DELTA).GT.THETN) THEN
            THETA=THETA+DELTA
	    CALL INTERP(RN,THETN,RN1,THETN1,THETA,R)
	    CALL IPOLAR(R,THETA,U,V)
	    WRITE(17,*) U,',',V
	    NINT=NINT+1
	    GOTO 80
	  END IF
	  RN=RN1
	  THETN=THETN1
	  GOTO 50
	ELSE
	  THETN=THETN1
	  RN=RN1
	  GOTO 60
	END IF
C
C
C	Interpolation on last interval and stop:
C
 600	THETN1=ALPHAN
 	RN1=RHON
C
 90	IF ((THETA.GT.THETN).AND.
     *    (THETA.LT.THETN1)) THEN
C
	IF ((THETN1-THETA).LT.DELTA) THEN
C
C***	  First check number of points written:
C
	  IF (NINT.NE.M) THEN
	    WRITE(9,91) INPUT NINT
	    WRITE(25,91) INPUT,NINT
 91	    FORMAT('# points in ',A17,' do not match with ',I4)
	  END IF
C*
	WRITE(9,*) '. . .Number of points written:',NINT,' . . .'
C*	  
C
C	  Now, stop this sample and reading new file:
	  WRITE(9,*) ' '
	  WRITE(9,*) ' '
	  WRITE(9,*) '. . .Reading next file. . .'
	  GOTO 5
	END IF
C
 	  CALL INTERP(RN,THETN,RN1,THETN1,THETA,R)
	  CALL IPOLAR(R,THETA,U,V)
	  WRITE(17,*) U,',',V
	  NINT=NINT+1
	  THETA=THETA+DELTA
	  GOTO 90
	END IF
C
	CLOSE(16)
	CLOSE(17)
C
C	Read in new file from list and repeat:
C
	WRITE(9,*) ' '
	WRITE(9,*) ' '
	WRITE(9,*) '. . .Reading next file. . .'
	GOTO 5
C
 999	WRITE(9,*) 'All files were treated'
 	CLOSE(25)
C
	WRITE(9,*) 'Enter a key to finish'
	PAUSE789
 789	CONTINUE
C
 	STOP
	END
	



C***************************************************************************
	SUBROUTINE POLAR(US,VS,RS,THETAS)
C***************************************************************************
C
C	Converts cartesian coordinates US,VS into polar coordinates
C	RS and THETA. Origin of local coordinate systems is at 0,0.
C	The calling unit gives US and VS as input. RS and THETAS are
C	returned to the calling unit. THETAS in Radians.
C
	DOUBLE PRECISION US,VS,RS,THETAS
	DOUBLE PRECISION ARG,PI
C
	PARAMETER(PI=3.14159265)
C
	RS=SQRT(US**2+VS**2)
C
C	First special cases:
C
	IF (US.GT.0.AND.VS.EQ.0) THEN
	  THETAS=0
	ELSE IF (US.EQ.0.AND.VS.GT.0) THEN
	  THETAS=PI/2
	ELSE IF (US.LT.0.AND.VS.EQ.0) THEN
	  THETAS=PI
	ELSE IF (US.EQ.0.AND.VS.LT.0) THEN
	  THETAS=3*PI/2
	END IF
C
C	For all other cases:
C
	IF (US.GT.0.AND.VS.GT.0) THEN
	  ARG=ABS(VS/US)
	  THETAS=ATAN(ARG)
	ELSE IF (US.LT.0.AND.VS.GT.0) THEN
	  ARG=ABS(US/VS)
	  THETAS=ATAN(ARG)+PI/2
	ELSE IF (US.LT.0.AND.VS.LT.0) THEN
	  ARG=ABS(VS/US)
	  THETAS=ATAN(ARG)+PI
	ELSE IF (US.GT.0.AND.VS.LT.0) THEN
	  ARG=ABS(US/VS)
	  THETAS=ATAN(ARG)+3*PI/2
	END IF
C
	RETURN
	END
	
	
	
	



C***************************************************************************
	SUBROUTINE IPOLAR(RHO,PHI,US,VS)
C***************************************************************************
C
C	Subroutine to convert polar coordinates RHO,PHI into
C	cartesian coordinates US,VS. US and VS are returned to the
C	calling unit. PHI is in Radians.
C
	DOUBLE PRECISION RHO,PHI,US,VS
C
	US=RHO*COS(PHI)
	VS=RHO*SIN(PHI)
C
	RETURN
	END






	
C*******************************************************************************	
	SUBROUTINE EXTREMA(XMIN,YMIN,XMAX,YMAX,I,INPUT)
C*******************************************************************************
C
C	Determines the minima and maxima of X,Y-pairs of values in
C	file INPUT. XMIN,XMAX,YMIN,YMAX are returned to the calling unit.
C	It also determines the number of points I of the input outline.
C
	INTEGER I
	DOUBLE PRECISION XMIN,YMIN,XMAX,YMAX
	CHARACTER*17 INPUT
C
	OPEN(15,FILE=INPUT,STATUS='OLD')
C
C	Determination of Xmin,Xmax,Ymin and Ymax:
C	XMIN, and counting the number of points (=I) in the input file:
C
	READ(15,*) XMIN,Y
	I=1
 5	READ(15,*,END=100) X,Y
 	I=I+1
 	IF (X.LT.XMIN) THEN
	  XMIN=X
	END IF
	GOTO 5
C

 100	REWIND(15)
C
C	XMAX:
C
	READ(15,*) XMAX,Y
 10	READ(15,*,END=200) X,Y
	IF (X.GT.XMAX) THEN
	  XMAX=X
	END IF
	GOTO 10
C
 200	REWIND(15)
C
C	YMIN:
C
	READ(15,*) X,YMIN
 15	READ(15,*,END=300) X,Y
	IF (Y.LT.YMIN) THEN
	  YMIN=Y
	END IF
	GOTO 15
C
 300	REWIND(15)
C
C	YMAX:
C
	READ(15,*) X,YMAX
 20	READ(15,*,END=400) X,Y
 	IF (Y.GT.YMAX) THEN
	  YMAX=Y
	END IF
	GOTO 20
C
 400	CLOSE(15)
C
 	RETURN
 	END



C*******************************************************************************
	SUBROUTINE CENTR(XCTR,YCTR,INPUT)
C*******************************************************************************
C
C	Subroutine to calculate the centroid (=point of gravity) of
C	a closed curve by trapezoidal approximation (after Davis, 1986,
C	pp. 346-347).
C
C	Input file is INPUT with the X- and Y coordinates of the points.
C	The subroutine returns the coordinates XCTR,YCTR of the centroid.
C
	DOUBLE PRECISION XN,YN,XN1,YN1,XCTR,YCTR
	DOUBLE PRECISION SUMA,SUMX,SUMY,AN,BN
	CHARACTER*17 INPUT
C
	OPEN(15,FILE=INPUT,STATUS='OLD')
C
C	Initialization and reading first point:
C
	SUMA=0.0
	SUMX=0.0
	SUMY=0.0
C
	READ(15,*) XN,YN
C
C	Memorize coordinates of starting point:
C
	AN=XN
	BN=YN
C
C	Read next point and calculate sums:
C
 100	READ(15,*,END=200) XN1,YN1
 	SUMA=SUMA+(YN1+YN)*(XN-XN1)/2
	SUMX=SUMX+(XN1**2+XN1*XN+XN**2)*(YN1-YN)/6
	SUMY=SUMY+(YN1**2+YN1*YN+YN**2)*(XN-XN1)/6
	XN=XN1
	YN=YN1
	GOTO 100
C
C	Calculate last interval (i.e.closing curve):
C
 200	XN1=AN
 	YN1=BN
C
 	SUMA=SUMA+(YN1+YN)*(XN-XN1)/2
	SUMX=SUMX+(XN1**2+XN1*XN+XN**2)*(YN1-YN)/6
	SUMY=SUMY+(YN1**2+YN1*YN+YN**2)*(XN-XN1)/6
C
	XCTR=SUMX/SUMA
	YCTR=SUMY/SUMA
C
	CLOSE(15)
C
	RETURN
	END





C***************************************************************************
	SUBROUTINE ROTATE(QN,OMEGN,QN1,OMEGN1,OMEGA,Q)
C***************************************************************************
C
C	Rotates OMEGN,OMEGN1 and OMEGA so, that OMEGN1 becomes zero
C	(i.e. all points are within the first quadrant of the cartesian
C	coord-system. It then determines the value of Q by linear 
C	interpolation, which is returned to the calling unit.
C	All angles are in Radians.
C
	DOUBLE PRECISION QN,OMEGN,QN1,OMEGN1,OMEGA,Q
	DOUBLE PRECISION ROT,PHIN,PHIN1,PHI
	DOUBLE PRECISION PI
C
	PARAMETER(PI=3.14159265)
C
	ROT=2*PI-OMEGN1
	PHIN=OMEGN+ROT
	PHIN1=0.0
	PHI=OMEGA-OMEGN1
C
C	Now interpolation between PHIN and PHIN1 to get the length of Q:
C
	CALL INTERP(QN,PHIN,QN1,PHIN1,PHI,Q)
C
	RETURN
	END





C***************************************************************************
	SUBROUTINE INTERP(RHON,PHIN,RHON1,PHIN1,PHI,RHO)
C***************************************************************************
C
C	Subroutine to linearly interpolate between two succeeding points.
C	Points enter the subroutine in polar coordinates. Interpolated
C	values are returned to the calling unit in polar coodinates.
C	All angular values are in Radians.
C
	DOUBLE PRECISION RHON,PHIN,RHON1,PHIN1,PHI,RHO
	DOUBLE PRECISION UN,VN,UN1,VN1,M
	DOUBLE PRECISION PI
C
	PARAMETER(PI=3.14159265)
C
	UN=RHON*COS(PHIN)
	VN=RHON*SIN(PHIN)
	UN1=RHON1*COS(PHIN1)
	VN1=RHON1*SIN(PHIN1)
	M=(VN1-VN)/(UN1-UN)
C
C	Interpolation of U and V (at angle PHI) from neighboring points
C	UN,VN and UN1,VN1:
C
C	First consider special cases for PHI:
C
	IF (PHI.EQ.O) THEN
	  U=(M*UN-VN)/M
	  V=0.0
	  GOTO 10
	ELSE IF (PHI.EQ.PI/2) THEN
	  U=0.0
	  V=-M*UN+VN
	  GOTO 10
	ELSE IF (PHI.EQ.PI) THEN
	  U=(-M*UN+VN)/M
	  V=0.0
	  GOTO 10
	ELSE IF (PHI.EQ.3*PI/2) THEN
	  U=0.0
	  V=M*UN-VN
	  GOTO 10
	ELSE IF (PHI.EQ.2*PI) THEN
	  U=(M*UN-VN)/M
	  V=0.0
	  GOTO 10
	END IF
C
C	All other cases:
C
	V=(VN-M*UN)*TAN(PHI)/(-M+TAN(PHI))
	U=V/TAN(PHI)
C
 10	RHO=SQRT(U**2+V**2)
C
	RETURN
	END
	
	
	
	
C***************************************************************************	
	SUBROUTINE DIAM(DMAX,ALPHA,BETA)
C***************************************************************************
C
C	Calculates the maximum diameter of an irregularly
C	shaped, closed plane curve.
C
C	Input file: POL_COORD, allready open from calling unit.
C	Coordinate system: Polar coordinates, with angular values
C	increasing in counterclock direction (=positive direction).
C
C	Output: The subroutine returns the maximum diameter (DMAX)
C	and the corresponding angular arguments ALPHA and BETA (in radians)
C	of the two opposite rays to the calling unit.
C
C	The programs first reads a a particular point PHI1,S1 of the
C	outline and then calculates the argument PHI2 of the oposite
C	point (=PHI2,S2) at 180¡ angular distance. Then, it reads
C	all succeeding points (PHI,S) of the outline, until PHI
C	exceeds the argument of PHI2. The sum of the rays D=S1+S2 is an
C	approximation of the diameter of the curve at this angular
C	value. The program calculates the maximum value DMAX of these
C	distances around the perimeter.
C
C	This approximation is good enough if the outline consists of a
C	large number of points (for example at 1000 points, the angular
C	distance between two subsequent points is 0.36¡, i.e. giving
C	an error for the angular estimate of the diameter of about ± 0.18¡).
C
C	NOTE: All angular values are in radians.
C
C
	INTEGER I,J
	DOUBLE PRECISION PI
	DOUBLE PRECISION D,DMAX
	DOUBLE PRECISION PHI,PHI1,PHI2,S,S1,S2
	DOUBLE PRECISION ALPHA,BETA
C
	PARAMETER(PI=3.14159265)
C
C	J is the number of the pair of polar coordinates in the input
C	file. J=1 is the first point, J=2 the second, etc.
C
C	Initialize DMAX and J
C
	DMAX=0.0
	J=1
C
C	Reading starting point of a particular diameter at argument PHI1:
C
 15	DO 20,I=1,J
 	  READ(16,*,END=999) PHI1,S1
 20	CONTINUE
C
C	Calculate argument at opposite position (=180¡ angular distance):
C
 	PHI2=MOD(PHI1+PI,2*PI)
C
C	Searching point to complete the diameter to 180¡ distance from
C	starting point. Then calculate diameter D and maximize.
C
 25	READ(16,*,END=999) PHI,S
C
C	Case I: Opposite rays do not cross the 360¡ ray (= positive x-axis):
C
 	IF (PHI2.GT.PHI1) THEN
	  IF (PHI-PHI2.GE.0.0) THEN
	    S2=S
	    D=S1+S2
	    IF (D.GT.DMAX) THEN
	      DMAX=D
	      ALPHA=PHI1
	      BETA=PHI
	    END IF
	    REWIND(16)
	    J=J+1
	    GOTO 15
	  END IF
	  GOTO 25
	  
	ELSE
C
C	Case II: Opposite rays do cross the 360¡ ray (= positive x-axis):
C
	  IF ((PHI.LE.PI).AND.((PHI-PHI2).GE.0.0)) THEN
	    S2=S
	    D=S1+S2
	    IF (D.GT.DMAX) THEN
	      DMAX=D
	      ALPHA=PHI1
	      BETA=PHI
	    END IF
	    REWIND(16)
	    J=J+1
	    GOTO 15
	  END IF
	  GOTO 25
	END IF
C
 999	RETURN
	END
