	program KeelWidth
C
C	Version 13 by Michael Knappertsbusch
C
C****	Important note:
C	KeelWidth13.out is for usage in combination with "Cycle II" programs Stat_Prep5.0.out,
C	Compose6.0.out and Stats9.0.out. Results from KeelWidth13.out can be plotted with XYPlot3.1.out.
C	Make sure you use the correct versions of these programs for your analysis !
C
C
C	First version from 25.11.2002.
C	Modified 29.11.2002: geological age added, and formatting slightly changed.
C	Modified  5.12.2002: geological age and various morphological parameters are read from
C	                     external file called MEASUREMENTS, which is obtained from Sprep53.out
C	Modified  6.12.2002: Formatting adapted.
C	Modified  9.12.2002: Calculation of keel angles added.
C	Modified 11.12.2002: Improvement of subroutine angle and subroutine correct added.
C	Modified 12.-13.12.2002: Corrections in subroutines EXTREMA and CORRECT.
C	Version 9:
C	Version 8 modified to version 9 (30.1.-28.2.2003): Addition of subroutines FPOINTS and AREA
C			     to calculate the "keel residual area" (=keel viw area minus areas of
C			     the four triangles (1) through (4).
C	Modified to Version 9.3, 19.12.2003: NCHAMB written to output as A2 instead of A1 character.
C	Modified to Version 10, 29.3.2005: Spiral angle is written to the output file and adaption
C			     to programs Stat_Prep30.out, Compose50.out, and Stats70.out.
C	Modified to Version 11.1 from version 10, 5.9.2011: Subroutine KREGIO implemented, so that 
C			     upper and lower leel region coordinates are written to two separate files.
C			     The coordinates of the upper keel region are written to file filename_UpKR,
C			     the coordinates of the lower keel region are written to file filename_LoKR.
C			     These newly generated files are then read by subroutine CIRCFIT to calculate
C			     centers and radii of osculating circles, that approximate the keel regions.
C			     In version 11.1 this calculation is done sing the separate program
C			     FitCircle1.3.2.out, where CIRCFIT is implemented.
C	Modified to Version 12.1 from version 11.1, 6.10.2011: In version 12.1 subroutine CIRCFIT is
C			     directly implemented in program KeelWidth, so that no external program is
C			     needed anly longer to determine osculating circles in the upper and lower
C			     keel region. The radii of the upper and lower keel regions are written
C			     to file RESKW... . An extra file List_KR_circles is generated, that contains
C			     for each outline the center coordinates and radii of the upper and lower
C			     osculating circle. This file List_KR_circles together with the traced
C			     outlines (_T files)can be read by program XYPLOT3.1 for plotting outlines
C			     with the circles.
C	Modified to Version 13 from Version 12.1, 23.11.2012: Addition of spiral convexity A and
C			     umbilical convexity B.
C			     Convexity A is the horizontal distance between the spiral apex (left) and
C			     the vertical midline (see drawings in CORRECT). Convexity B is the
C			     horizontal distance between the midline and the extreme right point on the
C			     umbilical side. 
C
C			             A < B means, that the outline is umbilico-convex.
C			             A = B means, that the outline is biconvex.
C			             A > B means, that the outline is spiro-convex.
C
C			     The same method to characterize the convexity of outlines in keel view is
C			     also applied in AMOR version higher than 3.17.
C
C			     In addition, subroutine CORRECT was modified in KeelWidth version 13, i.e. the
C			     coordinates of Point A2 are now explicitely calculated. (In version 12.1 and
C			     older versions, point A2 was simply set equal to point A under the assumption,
C			     that the starting point of the outline sits always on or is close to the 240
C			     pixels horizontal line. This was approximately true when upright positioning of
C			     the shells on the monitor was performed by hand, so that the apex becomes located
C			     on the Y=240 pixels line. In AMOR such precision for positioning is often not
C			     achieved, so that the starting point of the outline does no longer exactly
C			     coincide with the apex any more). Therefore, the new explicit determination of
C			     A2 slightly influences the coordinates of midpoint MA, and angles PHI1, PHI2 and
C			     PHI3.
C
C			
C
C
C
C
C
C	!!! Processing in batch mode. !!!
C	
C
C	Given the x,y coordinates of the silhouette of a globorotalia in keel view, this program
C	determines the width of the keel at the lower and upper end of the profile.
C	The lower keelwidth of the profile is determined at Q % of the vertical extension of the shell
C	in keel view, and the upper keelwidth is determined at (100-Q) % of the vertical extension of
C	the shell in keel view. The basis (=0% of the vertical extension of the shell in keel view) is
C	set at the lower keel.
C
C	The keelwidth is defined as the horizontal (x) distance between two points with identical
C	(or nearly) identical y-coordinates at Q % or (100-Q) % of the vertical (y) extension of
C	the shell in keel view.
C
C	The lower keelwidth (e.g. at Q% of the vertical extension) is called Y10, the upper keelwidth 
C	(at (100-Q) % of the vertical extension) is called Y90.
C
C	Because Y10 and Y90 are relative parameters with resepect to the shell size, it is NOT necessary
C	to use normalized outlines, so that the x,y coordinates can directly be taken from the traced
C	(_T) files.
C
C	The program calculates also the "keel-angles", and the "residual keel areas".
C
C
C	Input files:
C		Input is the file called "MEASUREMENTS", which contains the names of the specimens to
C		be treated. "MEASUREMENTS" must first be generated with program Sprep53.out.
C		From "MEASUREMENTS" the program reads the specimen names, which are the traced
C		(_T) outlines of the shells. As said above, it is not necessary to normalize the traced
C		files.
C		The other input files are the traced (_T) files containing the individual x,y coordinates
C		for each specimen.
C		Note, that the sequence of outline points needs to be in counterclockwise direction (in
C		program Trace outline extraction from the image occurred also in counterclockwise direction
C		with the starting point being on the y=240 pixels line or close to it.
C		The apex of the shell is always directing to the left hand side.
C		The _T files are read and treated sequentially.
C
C	Output files:
C		There is one main output file: From the name of the (_T) files, the species name, and the
C		number of chambers in the last whorl and the coiling direction (both entered by the
C		operator) the program calculates the corresponding name of the output file.
C
C		Example1:
C		Reading the specimen 502A106125K0001_T from "MEASUREMENTS" and getting menardA as species
C		code, 6chambers in the last whorl and a sinistral coiling direction, the output filename,
C		which contains the results will be RESKW_502A106125K_6Ch_s.
C
C		Two additional files filename_UpKR and filename_LoKR are written per outline containing the
C		cartesian coordinates of the upper and lower keel regions, respectively. The input filename
C		is used to compose the names of these two files.
C
C		Example2:
C		Reading the specimen 502A106125K0001_T from "MEASUREMENTS" the output file for the upper
C		keel region will be 502A106125K0001_UpKR.
C
C		One last file List_KR_circles is written containing the center coordinates and radii
C		of keel region osculating circles. Two lines are written for every outline. The first
C		line in the pair is for the lower keel region, the second line is for the upper keel region.
C
C		
C
C
C
C
C	List of subroutines:
C	EXTREMA, KWIDTH10, KWIDTH90, ANGLE, CORRECT, MIDPOINT, FPOINTS, ARSEG, KREGIO, CIRCFIT, CONVEX
C
C	List of functions:
C	TRIANGL
C
C
C
C
C
	INTEGER IERRL,IERRU			! Error handling: If CIRCFIT encounters a division by zero 
C						! during calculation of the Lower or Upper keel circle,
C						! IERRL and IERRU are set to 1, respectively.
	INTEGER N,NPOINTS,LAN			! LAN is a variable for the composition of the output files.
	INTEGER SEL				! SEL allows to select for output of Convexities A and B: if SEL=0, then
C						!   output of ARESLO,ARESUP; if SEL=1, then output of convexities A and B
C						!   instead of ARESLO and ARESUP.
	LOGICAL FILEX				! To check whether a file exists or not.
	REAL AGE,Q,R				! Age of the sample in Ma. R=¶X/¶Y. Q determines the relative position
C						! of D10 and D90 along the vertical axis of the keel view:
C						! For example, Q is 0.1 for D10 at 10% and D90 at 90% of the vertical extension. 
	REAL XCENTL,YCENTL,RL			! X and Y coordinates and radius for the Lower keel region circle.
	REAL XCENTU,YCENTU,RU			! X and Y coordinates and radius for the Upper keel region circle.
	DOUBLE PRECISION XMIN,YMIN,XMAX,YMAX	! XMIN,XMAX are the leftmost or rightmost X-coordinates of the outline,
C						! YMIN,YMAX are the lowermost or uppermost Y-coordinates of the outline.
	DOUBLE PRECISION YA,XB,YC,XD		! X, and Y-coordinates of the corresponding values for the extrema (see
C						! in subroutine EXTREMA).
	DOUBLE PRECISION Y10,Y90,D10,D90	! Y10,Y90 are the y-coordinates at the lower or upper keelwidths, D10 
	DOUBLE PRECISION D9010,DMIN,DMAX	! and D90 are the values of the lower and upper keelwidths, D9010 is
	DOUBLE PRECISION DX,DY			! the absolute difference between D90 and D10, DMIN is the smaller and DMAX
	DOUBLE PRECISION XL10,XR10,XL90,XR90	! DMAX is the larger of D10 and D90, and DX is the difference between
C						! XMAX and XMIN (= the spiral height of the shell).
	DOUBLE PRECISION DELTAX,DELTAY,AREA	! DELTAX (=¶X), DELTAY (=¶Y): DELTAX is identical to DX, and DELTAY
C						! is identical to DY; DX and DY are values calculated values, while
C						! DELTAX and DELTAY are taken from file MEASUREMENTS.
	DOUBLE PRECISION AREA1,AREA2		! Areas (1) to (4) necessary for determination of "keel residual areas".
	DOUBLE PRECISION AREA3,AREA4		
	DOUBLE PRECISION XSI1,XSI2,XSI3,XSI4	! Areas of triangles (1) to (4) for determination of "keel residual areas". 
	DOUBLE PRECISION ARESLO,ARESUP		! Lower and upper "keel residual areas".
	DOUBLE PRECISION XMA,YMA,XMB,YMB	! XMA,YMA: Midpoint coords between A and A1; XMB,YMB: Midpoint between B and B1;
	DOUBLE PRECISION XMC,YMC,XMD,YMD	! XMC,YMC: Midpoint coords between C and C1; XMD,YMD: Midpoint between D and D1.
	DOUBLE PRECISION XF1,YF1,XF2,YF2	! Coords of the two "F-points" F1 and F2.
	DOUBLE PRECISION PHI1,PHI2		! Upper and lower keel angle (PHI1 and PHI2, respectively), in degrees.
	DOUBLE PRECISION PHI3			! PHI3 is the spiral angle in keel view, in degrees.
	DOUBLE PRECISION TRIANGL		! Function TRIANGL, that determines the area of the triangle given by 3 points.
	DOUBLE PRECISION DEV,SUM1		! DEV: Control variable (in %) for estimation of deviation of sum of (AREA1+AREA2+AREA3+AREA4)
C						! from AREA (given in file "MEASUREMENTS" during determination of ARESLO and ARESUP.
	DOUBLE PRECISION A,B,RATIO		! A indicates the spiral convexity a, b indicates the umbilical convexity b,
C						! and RATIO is A/B.
C
	CHARACTER*66 HEADER			! Header line in file MEASUREMENTS
	CHARACTER*7 SPECIES			! Species name (same name encoding as used in programs KeelBend, Homolo, etc).
	CHARACTER*17 INPUT			! Name of the file with the outline coordinates (_T files, 17 chars long).
	CHARACTER*31 OUTPUT			! Name of the file with the results.
	CHARACTER*1 COIL			! Coiling direction.
	CHARACTER*2 NCHAMB			! NCHAMB: Number of chambers.
	CHARACTER*20 NUPKR,NLOKR		! Names of the files containing the X,Y coords of the upper (NUPKR) and the lower (NLOKR) keel regions
C
C
C
C
	WRITE(9,*) '***************************************************'
	WRITE(9,*) '*                                                 *'
	WRITE(9,*) '*             Program KeelWidth                   *'
	WRITE(9,*) '*  determines width of the Keel in silhouettes    *'
	WRITE(9,*) '*  of globorotaliid foraminifers in side view     *'
	WRITE(9,*) '*                                                 *'
	WRITE(9,*) '*         !!! Batch processing mode !!!           *'
	WRITE(9,*) '*                                                 *'
	WRITE(9,*) '*     Version 13  By Michael Knappertsbusch       *'
	WRITE(9,*) '*                                                 *'
	WRITE(9,*) '*                                                 *'
	WRITE(9,*) '*     For usage in combination with programs      *'
	WRITE(9,*) '*      Stat_Prep5.0.out, Compose6.0.out and       *'
	WRITE(9,*) '*                   Stats9.0.out                  *'	
	WRITE(9,*) '*                                                 *'
	WRITE(9,*) '***************************************************'
	WRITE(9,*) ' '
	WRITE(9,*) ' '
C
C
C
C
C	Recording errors during program run:
C
	OPEN(25,FILE='Errors_during_KeelWidth',STATUS='NEW')
C
	WRITE(9,*) '. . .Enter KeelWidth (Q). . .'
	WRITE(9,*) ' (Q is between 0 and 1)'
	WRITE(9,*) 'Recommended Q value for menardiforms: 0.05)'
	READ(9,*) Q
	WRITE(9,*) ' '
	WRITE(9,15) Q*100.0,100*(1.0-Q)
 15	FORMAT('Lower (=D10) and upper KeelWidth (=D90)
     1 at',F5.1,'% and ',F5.1,'% of axial extension')
C
	WRITE(9,*) ' '
	WRITE(9,*) '. . .Output of ARESLO,ARESUP or 
     1 output of convexities A,B ?. . .'
	WRITE(9,*) '   enter 0 or 1'
	READ(9,*) SEL
C
C
C	Get filenames from file list with name "MEASUREMENTS", open that file and
C	read the names of the files with the traced outlines.
C	Note: The file MEASUREMENTS must first be generated with program Sprep53.out
C
	OPEN(20,FILE='MEASUREMENTS',STATUS='OLD')
C	Determine sample name from INPUT, and write the sample name into the
C	header of file Errors_during_KeelWidth:
C
	READ(20,'(F5.2)') AGE
	READ(20,'(A66)') HEADER
	READ(20,6) INPUT,AGE,NPOINTS,DELTAX,DELTAY,AREA
 6	FORMAT(A17,X,F5.2,X,I4,X,F6.1,X,F6.1,X,F7.4)
C
	WRITE(25,3) INPUT(1:11)				! 25 ='Errors_during_KeelWidth'
 3	FORMAT('Sample name: ',A11)
	REWIND(20)
C
C***	The following 2 lines are necessary to reposition the cursor
C	to the beginning of the first data record in file MEASUREMENTS:
C
	READ(20,'(F5.2)') AGE
	READ(20,'(A66)') HEADER
C
C
C****	Preparations for output of results:
C	Write header line for results to screen:
C
     	WRITE(9,1)
 1	FORMAT('The values of D10, D90, D9010,
     1 and DMIN are written as relative values
     2 with respect to DX (e.g. the spiral height
     3 of the shell).')
C
C
C***	Get species names, # of chambers and the coiling direction from screen:
C
	WRITE(9,*) 'Enter species name (7chars),
     1 # of chambers (2 chars), and the coiling
     2 (s or d or n [for none] ), separated by commas'	! choose n if the coiling is not specified
     	WRITE(9,*) 'If # of chambers or if coiling is
     1 not specified enter n (for none)
     2 for each variable'
     	WRITE(9,*) 'Example: menardA, 6,s'
	READ(9,14) SPECIES,NCHAMB,COIL
 14	FORMAT(A7,1X,A2,1X,A1)
C
	WRITE(9,*) ' '
	IF (SEL.EQ.0) THEN
	  WRITE(9,12) SEL
	  WRITE(9,2)				! Labelling columns to screen if SEL=0 (output with ARESLO, ARESUP, and DEV)
	END IF
C
	IF (SEL.EQ.1) THEN
	  WRITE(9,13) SEL
	  WRITE(9,7)				! Labelling columns to screen if SEL=1 (output with convexities A, B, and RATIO)
	END IF
C
C
 2	FORMAT('Specimen,            Age
     1,#pts,Species,#Ch,Coil,¶X (µm),¶Y (µm)
     2,R=¶X/¶Y,Ar (mm2),  D10 %,   D90 %,   D9010%
     3,  Dmin%,   Dmax%,   Phi1¡,   Phi2¡
     4,  ARESLO (mm2),  ARESUP (mm2), DEV (in %)
     5, Phi3¡,  RL (µm),  RU (µm)')
C
 7	FORMAT('Specimen,            Age
     1,#pts,Species,#Ch,Coil,¶X (µm),¶Y (µm)
     2,R=¶X/¶Y,Ar (mm2),  D10 %,   D90 %,   D9010%
     3,  Dmin%,   Dmax%,   Phi1¡,   Phi2¡
     4,  A (µm),        B (µm),    Ratio A/B
     5,     Phi3¡,  RL (µm),  RU (µm)')
C
 12	FORMAT(I1,' (Output with ARESLO, ARESUP, DEV)')
 13	FORMAT(I1,' (Output with convexities A, B, and RATIO)')
C
C	Now, compose the name of the output file:
C
	LAN=2						! Sets the value for NCHAMB to 2 characters
	IF (NCHAMB.EQ.'n') THEN
	OUTPUT='RESKW_'//INPUT(1:11)//'_'
     1//SPECIES//'_nCH_'//COIL
     	  LAN=1
	ELSE IF (NCHAMB.EQ.' 1') THEN
	OUTPUT='RESKW_'//INPUT(1:11)//'_'
     1//SPECIES//'_1CH_'//COIL
     	  LAN=1
	ELSE IF (NCHAMB.EQ.' 2') THEN
	OUTPUT='RESKW_'//INPUT(1:11)//'_'
     1//SPECIES//'_2CH_'//COIL
     	  LAN=1
	ELSE IF (NCHAMB.EQ.' 3') THEN
	OUTPUT='RESKW_'//INPUT(1:11)//'_'
     1//SPECIES//'_3CH_'//COIL
     	  LAN=1
	ELSE IF (NCHAMB.EQ.' 4') THEN
	OUTPUT='RESKW_'//INPUT(1:11)//'_'
     1//SPECIES//'_4CH_'//COIL
     	  LAN=1
	ELSE IF (NCHAMB.EQ.' 5') THEN
	OUTPUT='RESKW_'//INPUT(1:11)//'_'
     1//SPECIES//'_5CH_'//COIL
     	  LAN=1
	ELSE IF (NCHAMB.EQ.' 6') THEN
	OUTPUT='RESKW_'//INPUT(1:11)//'_'
     1//SPECIES//'_6CH_'//COIL
     	  LAN=1
	ELSE IF (NCHAMB.EQ.' 7') THEN
	OUTPUT='RESKW_'//INPUT(1:11)//'_'
     1//SPECIES//'_7CH_'//COIL
     	  LAN=1
	ELSE IF (NCHAMB.EQ.' 8') THEN
	OUTPUT='RESKW_'//INPUT(1:11)//'_'
     1//SPECIES//'_8CH_'//COIL
     	  LAN=1
	ELSE IF (NCHAMB.EQ.' 9') THEN
	OUTPUT='RESKW_'//INPUT(1:11)//'_'
     1//SPECIES//'_9CH_'//COIL
     	  LAN=1
	END IF
C
	IF (LAN.GT.1) THEN
	OUTPUT='RESKW_'//INPUT(1:11)//'_'
     1//SPECIES//NCHAMB//'CH_'//COIL
     	END IF
C
C
	OPEN(18,FILE='List_KR_circles',STATUS='NEW')	! File List_KR_circles contains the X,Y center coordinates and
C							! radii of the Lower and Upper osculating keel region circles.
C							! For each outline there is a pair of lines: The first line
C							! holds the values for the Upper circle, the second line is
C							! for the Lower circle.
	WRITE(18,4)					! Write header line in file List_KR_circles
 4	FORMAT('File with keelregion  Xcenter
     1  Ycenter   Radius')
C
	OPEN(35,FILE=OUTPUT,STATUS='NEW')		! File RESKW... contains the results
	WRITE(35,*) 'Notes:'				! Write various notes to RESKW...
	WRITE(35,15) Q*100.0,100*(1.0-Q)
	WRITE(35,1)					 
	WRITE(35,*) ' '
	IF (SEL.EQ.0) THEN
	  WRITE(35,12) SEL				! A "0" informs programs StatPrep, Compose, and Stats which format to use
	  WRITE(35,2)					! Labelling columns in file RESKW... if SEL=0 (output with ARESLO, ARESUP, DEV)
	END IF
C
	IF (SEL.EQ.1) THEN
	  WRITE(35,13) SEL				! A "1" informs programs StatPrep, Compose, and Stats which format to use
	  WRITE(35,7)					! Labelling columns in file RESKW... if SEL=1 (output with convexities A, B, and RATIO)
	END IF
C
C
C
C
C
C
C******	Main loop starts at statement number 5:
C
 5	READ(20,6,END=999) INPUT,AGE,NPOINTS,DELTAX,DELTAY,AREA	! 	20 = 'MEASUREMENTS'
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				! 25 = 'Errors_during_KeelWidth'
 8	  FORMAT('File ',A17,' is missing')
	  WRITE(9,*) '. . .Program continues. . .'
 	  GOTO 5
	END IF
C
C
C
C	
C*****	Calculations for each INPUT file start here:
C	First, determination of extrema of outline coordinates:
C
     	CALL EXTREMA(INPUT,XMIN,YA,XB,
     1 YMIN,XMAX,YC,XD,YMAX,N)
C
C-----------------------------------------------CHECK: EXTREMA returns correct coordinates at this point (27.11.2012).
C
	DX=XMAX-XMIN				! DX is the maximum extension of the shell in X direction (=spiral height ¶X).
	DY=YMAX-YMIN				! DY is the maximum extension of the shell in Y direction (=axial diameter ¶Y).
	R=DX/DY
C
C*****	Determination of the y-coordinates of the lower (Y10) and upper (Y90) keelwidths:
	Y10=YMIN+Q*(YMAX-YMIN)			! Y-position of lower keelwidth (recommended: Q=0.05)
	Y90=YMIN+(1.0-Q)*(YMAX-YMIN)		! Y-position of upper keelwidth
C
C
C*****	Added 31.8.2011:
C	Write upper and lower keel region coordinates to two separate external files. These files are then
C	sent to subroutine CIRCFIT to calculate osculating circles that describe the
C	keel regions.
C
	CALL KREGIO(Y10,Y90,INPUT,NLOKR,NUPKR)	! Returns the file names for the lower and upper keel regions
	CALL CIRCFIT(NLOKR,XCENTL,YCENTL,RL,IERRL)	! Returns Lower keel region circle
	CALL CIRCFIT(NUPKR,XCENTU,YCENTU,RU,IERRU)	! Returns Upper keel region circle
C
C***	Write circles to external file List_KR_circles. For each outline a pair of lines are written: The first
C	line contains the Lower keel region circle, the second line contains the Upper keel region circle.
C	Together with the traced (_T) outlines the file List_KR_circles can then be read by program XYPlot3.1 to
C	plot outlines with the keel region circles on the computer screen.
C
	IF (IERRL.EQ.1) THEN
	  WRITE(18,21) INPUT(1:15)//'_LoKR'
	ELSE
	  WRITE(18,22) INPUT(1:15)//'_LoKR',XCENTL,YCENTL,RL
	END IF
C
	IF (IERRU.EQ.1) THEN
	  WRITE(18,21) INPUT(1:15)//'_UpKR'
	ELSE
	  WRITE(18,22) INPUT(1:15)//'_UpKR',XCENTU,YCENTU,RU
	END IF
C
 21	FORMAT(A20,' Error: Division by zero in CIRCFIT.
     1 Modify input coordinates and retry for this file')
 22	FORMAT(A20,3(1X,F8.3))
C
C
C
C*****	Now determine corresponding X coordinates on left and right limbs of shell silhouette for
C	upper and lower keelwidths:
C
	CALL KWIDTH10(Y10,XL10,XR10,D10,INPUT)	! Lower keelwidth (at RATIO (default 10)% of height of profile)
	CALL KWIDTH90(Y90,XL90,XR90,D90,INPUT)	! Upper keelwidth (at RATIO (default 90)% of height of profile)
C	
C****	The so determined lower and upper keelwidths are indicative, whether the test is quasi-symmetrical
C	in axial view. However, if the test is slightly "twisted", as is often the case for example
C	in Globorotalia tumida specimens, only the minimum value of the two parameters is indicative,
C	with the larger value being the result of a non-ideal axial view. Therefore, the minimum
C	value DMIN is determined and written to the result file. 
C
	DMIN=MIN(D10,D90)
	DMAX=MAX(D10,D90)
C
C
C	Calculate difference between D90 and D10 (indcative for the degree of asymmetry of
C	the shell in axial view.
C
	D9010=ABS(D90-D10)
C
C
C****	Calculation of the upper (PHI1) and lower (PHI2) keel angles:
C	First determine correct vectors to midpoints, that determine PHI1 and PHI2 with
C	subroutine CORRECT:
C
C
	CALL CORRECT(INPUT,XMIN,YA,XB,
     1 YMIN,XMAX,YC,XD,YMAX,XMA,YMA,XMB,YMB,
     2 XMC,YMC,XMD,YMD)
C
C-----------------------------------------------CHECK: CORRECT returns correct coordinates at this point (27.11.2012).
C
C
C	"Midpoints" are now points MA, MB, MC, and MD with the coordinates (XMA,YMA),
C	MB(XMB,YMB), MC(XMC,YMC), and MD(XMD,YMD).
C
C	Now, calculate spiral convexity A and umbilical convexity B from midpoints MA,
C	MB, MC, and MD (added for KeelWidth13 on 23.11.2012):
C
C
	CALL CONVEX(XMA,XMB,XMC,XMD,A,B)			! Returns A and B
	RATIO=A/B
C
C-----------------------------------------------CHECK: CONVEX returns correct values for A and B at this point (28.11.2012).
C
C
C	Next, calculate PHI1, PHI2, AND PHI3 from midpoint vectors:
	
	CALL ANGLE(XMC,YMC,XMD,YMD,XMA,YMA,PHI1)		! Returns PHI1
	CALL ANGLE(XMA,YMA,XMB,YMB,XMC,YMC,PHI2)		! Returns PHI2
	CALL ANGLE(XMD,YMD,XMA,YMA,XMB,YMB,PHI3)		! Returns PHI3
C
C
C
C****   Calculation of "keel residual area" (added 30.1.2003):
C
C	Note, that residual areas ARESLO and ARESUP are only written to file if SEL is set to 0.
C	If SEL=1, spiral convexities A and B are written to file instead.
C
C	First determination of coordinates of points F1 and F2 with
C	subroutine FPOINTS.
C
	CALL FPOINTS(XMA,YMA,XMB,YMB,XMC,YMC,
     1 XMD,YMD,XF1,YF1,XF2,YF2)
C
C	Now calculate the four areal segments (1) to (4):
	CALL ARSEG(1,XF1,YF1,INPUT,XMA,YMA,
     1  XMB,YMB,AREA1)						! Returns AREA1
	CALL ARSEG(2,XF2,YF2,INPUT,XMB,YMB,
     2  XMC,YMC,AREA2)						! Returns AREA2
	CALL ARSEG(3,XF2,YF2,INPUT,XMC,YMC,
     3  XMD,YMD,AREA3)						! Returns AREA3
	CALL ARSEG(4,XF1,YF1,INPUT,XMD,YMD,
     4  XMA,YMA,AREA4)						! Returns AREA4
C
C
C+++	Check if AREA1 through AREA4 are correctly calculated (in
C	percent of the area given in file "MEASUREMENTS":
	SUM1=(AREA1+AREA2+AREA3+AREA4)/1.E6
	DEV=100.*(SUM1-AREA)/AREA
C+++	End Check
C     
C
C	Next, calculate the areas PHI1 through PHI 4 of the four triangles
C	defined by points MA,MB,MC,MD,F1 and F2:
C
	XSI1=TRIANGL(XMA,YMA,XF1,YF1,XMB,YMB)			! Triangle (MA,F1,MB)
	XSI2=TRIANGL(XMB,YMB,XMC,YMC,XF2,YF2)			! Triangle (F2,MB,MC)
	XSI3=TRIANGL(XF2,YF2,XMC,YMC,XMD,YMD)			! Triangle (MD,F2,MC)
	XSI4=TRIANGL(XMA,YMA,XF1,YF1,XMD,YMD)			! Triangle (MA,F1,MD)
C
C	Finally, calculation of lower and upper "keel residual areas" in mm2:
C
	ARESLO=(AREA1-XSI1 + AREA2-XSI2)/1.E6			! Lower "keel residual area", in mm2
	ARESUP=(AREA3-XSI3 + AREA4-XSI4)/1.E6			! Upper "keel residual area", in mm2
C
C
C
C
C*****  Output of results to screen and to file RESKW....
C	In order to make D10, D90, D9010, DMIN and DMAX independent of size these values
C	are written as percentages of the spiral height DX of the shell.
C
C
	IF (SEL.EQ.0) THEN					! Output if option 0 is selected (output
C								  of ARESLO,ARESUP to file)
C
	WRITE(9,10) INPUT,AGE,NPOINTS,SPECIES,
     1 NCHAMB,COIL,DX,DY,R,AREA,D10/DX,D90/DX,
     2 D9010/DX,DMIN/DX,DMAX/DX,PHI1,PHI2,
     3 ARESLO,ARESUP,DEV,PHI3,RL,RU

	WRITE(35,10) INPUT,AGE,NPOINTS,SPECIES,
     1 NCHAMB,COIL,DX,DY,R,AREA,D10/DX,D90/DX,
     2 D9010/DX,DMIN/DX,DMAX/DX,PHI1,PHI2,
     3 ARESLO,ARESUP,DEV,PHI3,RL,RU
C
	END IF
C
C
	IF (SEL.EQ.1) THEN					! Output if option 1 is selected (convexities
C								  A, B and the ratio A/B are written to file
C								  in place of ARESLO and ARESUP)
C
	WRITE(9,30) INPUT,AGE,NPOINTS,SPECIES,
     1 NCHAMB,COIL,DX,DY,R,AREA,D10/DX,D90/DX,
     2 D9010/DX,DMIN/DX,DMAX/DX,PHI1,PHI2,
     3 A,B,RATIO,PHI3,RL,RU
       
	WRITE(35,30) INPUT,AGE,NPOINTS,SPECIES,
     1 NCHAMB,COIL,DX,DY,R,AREA,D10/DX,D90/DX,
     2 D9010/DX,DMIN/DX,DMAX/DX,PHI1,PHI2,
     3 A,B,RATIO,PHI3,RL,RU
C
	END IF
C
C
C
 10	FORMAT(A17,', ',F5.2,', ',I4,', ',A7
     1',',A2,', ',A1,', ',F6.1,', ',F6.1
     2', ',F6.4,', ',F7.4,', ',F7.2,', ',F7.2
     3', ',F7.2,', ',F7.2,', ',F7.2,', ',F7.3
     4', ',F7.3,', ',F7.4,',      ',F7.4,
     5',     ',F7.3,',     ',F7.3,', ',F8.3
     6', ',F8.3)
C
 30	FORMAT(A17,', ',F5.2,', ',I4,', ',A7
     1',',A2,', ',A1,', ',F6.1,', ',F6.1
     2', ',F6.4,', ',F7.4,', ',F7.2,', ',F7.2
     3', ',F7.2,', ',F7.2,', ',F7.2,', ',F7.3
     4', ',F7.3,', ',F7.2,',      ',F7.2,
     5',     ',F7.3,',     ',F7.3,', ',F8.3
     6', ',F8.3)
C
C
C****	Reading next file:
C
	GOTO 5
C
 999	WRITE(9,*) ' '
 	WRITE(9,*) 'All files treated'
C
C	Added 31.8.2011:
	WRITE(9,*) 'Separate files for
     1 upper and lower keel regions
     2 written'
     	  WRITE(9,*) 'One extra file "List_KR_circles"
     1 written for input to program XYPlot3.1.out'
C
 	PAUSE1000
 1000	CONTINUE
C
	STOP
	END
	








C*******************************************************************************
	SUBROUTINE CIRCFIT(FNAME,XC,YC,R,ERR)
C*******************************************************************************
C
C	Version 1.2 by Michael Knappertsbusch, 24.8.2011
C
C	Given are cartesian x,y coordinates of a finite series of points, for which
C	a circle, that best passes through the points is to be found.
C
C	(Comment added on 7.10.2011):
C	Input to the subroutine is a filename FNAME with the X,Y coordinates for which
C	a fitting circle is to be found. The center coordinates XC and YC, the radius R,
C	and an error variable ERR are returned to the calling unit .
C
C	The method applies a least-squares circle fit method described by R. Bullock.
C	URL: http://www.dtcenter.org/met/users/doc/write_ups/circle_fit.pdf
C
C	The input points are stored in the external file with the name INPUT
C	X,Y are the cartesian coordinates of the input points.
C
C
C	N     = Running number of points
C	NPTS  = Number of input points
C	X,Y   = Cartesian coordinates of input points
C	SX    = Sum of the X values
C	SY    = Sum of the Y values
C	MEANX = Mean of the X values
C	MEANY = Mean of the Y values
C	U,V   = Center-transformed values of X and Y, respectively
C	SU    = Sum of the U values
C	SV    = Sum of the V values
C	SUU   = SUM of U*U
C	SVV   = SUM of V*V
C	SUV   = Sum of U*V
C	SUUU  = Sum of U*U*U
C	SVVV  = Sum of V*V*V
C	SUVV  = sum of U*V*V
C	SVUU  = Sum of V*U*U
C	UC,VC = Coordinates of the center of the best-fit circle in the U,V coord'system
C	XC,YC = Coordinates of the center of the best-fit circle in the X,Y coord'system
C	R     = Radius of the best-fit circle
C	ERR   = Variable for error handling: There is an error if ERR=1, else ERR=0
C	
C
	INTEGER N,NPTS,ERR
	REAL X,Y,SX,SY,MEANX,MEANY
	REAL U,V,SU,SV,SUU,SVV,SUV
	REAL SUUU,SVVV,SUVV,SVUU
	REAL UC,VC,XC,YC,ALPHA,R
	REAL TERM1,TERM2,TERM3,TERM4,TERM5
	CHARACTER*20 FNAME
C
	ERR=0			! Initialize ERR to 0 (no error)
C
C
C
C***	First open the input file, read the X,Y coordinates and determine
C	the number NPTS of input points.
C
	OPEN(15,FILE=FNAME,STATUS='OLD')
C
C	Count number of points:
C
	N=0
 10	READ(15,*,END=20) X
 	N=N+1
	GOTO 10
C
 20	NPTS=N
 	REWIND(15)
C
C
C***	Calculate means MEANX and MEANY of X and Y, respectively:
C
	SX=0.0				! Initialize Sum SX
	SY=0.0				! Initialize Sum SY
 25	READ(15,*,END=30) X,Y
	SX=SX+X				! Summing up X values
	SY=SY+Y				! Summing up Y values
	GOTO 25
C
 30	MEANX=SX/NPTS
 	MEANY=SY/NPTS
C
	REWIND(15)
C
C
C***	Transform X and Y values to new coordinates U and V, respectively,
C	and write them to a temporary scratch file (file will no longer exist
C	after termination of program:
C
	OPEN(16,STATUS='SCRATCH')
 40	READ(15,*,END=50) X,Y
 	U=X-MEANX
	V=Y-MEANY
	WRITE(16,*) U,V
	GOTO 40
 50	CLOSE(15)
 	REWIND(16)
C
C
C***	Determine sums and sums of products of U and V values.
C	First initialize every summatiuon:
C
	SU=0.0
	SV=0.0
	SUU=0.0
	SVV=0.0
	SUV=0.0
	SUUU=0.0
	SVVV=0.0
	SUVV=0.0
	SVUU=0.0
C
 60	READ(16,*,END=70) U,V
	SUU=SUU+U*U
	SVV=SVV+V*V
	SUV=SUV+U*V
	SUUU=SUUU+U*U*U
	SVVV=SVVV+V*V*V
	SUVV=SUVV+U*V*V
	SVUU=SVUU+V*U*U
	GOTO 60
C
 70	CLOSE(16)		! Closes and deletes temporary scratch file
C
C
C***	Now calculate coordinates of the midpoint (UC,VC) of the fitting circle
C	(in the transformed U,V coordinate system)
C
	TERM3=2*(SUU*SUV*SUV-SVV*SUU*SUU)
	TERM5=2*(SUV*SUV-SVV*SUU)
C
C***	Error handling: Check for divisions by zero:
C
	IF ((SUU.EQ.0.0).OR.(TERM3.EQ.0.0).OR.
     1  (TERM5.EQ.0.0)) THEN
	  ERR=1
	  RETURN
	END IF
C
	TERM1=0.5*(SUUU+SUVV)/SUU
	TERM2=SUV*SUV*(SUUU+SUVV)-SUU*SUV*(SVVV+SVUU)		
	UC=TERM1-TERM2/TERM3
C
	TERM4=SUV*(SUUU+SUVV)-SUU*(SVVV+SVUU)
	VC=TERM4/TERM5
C
C
C***	Back-transform the center coordinates from the U,V coordinate
C	system into the original X,Y coordinate system:
C
	XC=UC+MEANX
	YC=VC+MEANY
C
C
C***	Determine the radius R of the best fit circle:
C
	ALPHA=UC**2+VC**2+(SUU+SVV)/NPTS
	R=SQRT(ALPHA)
C
C
	RETURN
 	END

















C*******************************************************************************
	SUBROUTINE KREGIO(V10,V90,FNAME,FLOKR,FUPKR)
C*******************************************************************************
C
C	By Michael Knappertsbusch, 31.8.2011
C
C	Given are the cartesian coordinates (U,V) of points of a plane closed curve
C	(globorotalid outline in profile view) and the V-coordinates of the lower
C	(V19) and upper (V90) keel widths.
C	This subroutine writes all U,V coordinates of the lower keelregion (V less or equal V10)
C	to file filename_LoKR. Similarly, it exports all coordinates of the upper
C	keel region (V greater or equal V90) to file filename_UpKR.
C	The two new files can then be imported to the separate program FitCircle1.3.out
C	for the determination of osculating circles (center coordinates, radii), that
C	describe the lower and upper keelregions.
C
C	Input to the subroutine are V10,V90 and FNAME.
C	The names of character variables FLOKR and FUPKR are returned to the calling unit.
C
C
	DOUBLE PRECISION V10,V90,U,V
	CHARACTER*17 FNAME
	CHARACTER*20 FLOKR,FUPKR		! Variables for the new file names
C
C*****	Determination of the output file names:
C
	FLOKR=FNAME(1:15)//'_LoKR'		! FLOKR contains the Filename for the LOwer KeelRegion
	FUPKR=FNAME(1:15)//'_UpKR'		! FUPKR contains the Filename for the UPper KeelRegion
C
	OPEN(30,FILE=FLOKR,STATUS='NEW')
	OPEN(31,FILE=FUPKR,STATUS='NEW')
	OPEN(32,FILE=FNAME,STATUS='OLD')
C
 10	READ(32,*,END=30) U,V
	IF (V.LE.V10) THEN			! Coordinates of lower keel region
	  WRITE(30,20) U,V
	ELSE IF (V.GE.V90) THEN			! Coordinates of upper keel region
	  WRITE(31,20) U,V
	END IF
 20	FORMAT(F8.3,',',F8.3)
	GOTO 10
C
 30	CLOSE(30)
 	CLOSE(31)
	CLOSE(32)
C
	RETURN
	END









C*******************************************************************************
	SUBROUTINE ARSEG(ID,XF,YF,INPUT,
     1 XSTART,YSTART,XEND,YEND,AREA)
C*******************************************************************************
C
C	By Michael Knappertsbusch, 27.1.2003
C	Latest modifications: 18.2.2003
C
C	Given the cartesian x,y coordinates of a plane closed curve, this
C	subroutine calculates the area of one of the four areal segments labelled with
C	(1), (2), (3), or (4). These areas are necessary to determine the "keel residual
C	areas" in the upper and lower keel region of a globorotalid foraminifer
C	when viewed in keel position (see figure below). AREA is returned to the calling
C	unit in units of input coordinates (mostly in micrometers).
C
C
C
C	   D1(XD1,YMAX)        MD  D(XD,YMAX)
C                          *....*...*
C                    .....      |   .
C                  ..     (4)   |    .   Outline of a keeled foram in keel position
C   A2(XMIN,YA2)  *             |     .
C   A(XMIN,YA) -->* (=240 PIXEL LINE)  .
C                 .             |       .
C           MA    *-------------+ F1    * C1(XMAX,YC1)
C                 .             |       .
C                 .             |  (3)  .
C    A1(XMIN,YA1) *    (1)      |       .
C                  .         F2 +-------* MC
C                   .           |       .
C                    .          |       .
C                     .         |       .
C                      .        |  (2)  * C(XMAX,YC)
C                       .       |      .
C                        :      |      .
C                         .     |     .
C                          .    |    .
C               B(XB,YMIN)  *...*...* B1(XB1,YMIN)
C                               MB
C
C
C
C	Input into the subroutine are the coordinates of one of the "F-points" (e.g.
C	F1 or F2), the name of the file that contains the x,y coordinates of the
C	outline (INPUT), the start and end points determining the particular segment
C	(example: for segment (1) the starting point is MA, the end point is MB and
C	the "F-point" is F1). The result (in units as given in the input file name)
C	is returned to the calling unit by the variable AREA.
C
C	The method of calculation of the area is by summing up the vector products
C	determined by two neighbouring points of the outline (e.g. Ai=0.5*ABS(ri X Ri+1) ).
C
C	Description of variables:
C		XF,YF = Coordinates of one of the "F-points"
C		XSTART,YSTART,XEND,YEND = Coordinates of start and end points of outline
C		                          segment, that determines one of the four
C					  areas (1) to (4).
C		AREA = Area of one of the four segments (1) through (4).
C		INPUT = File name with the x,y coordinates of the outline (20 chars).
C
C	List of Functions required:
C	TRIANGL
C
C
	INTEGER ID
	DOUBLE PRECISION XF,YF,XSTART,YSTART
	DOUBLE PRECISION XEND,YEND,AREA
	DOUBLE PRECISION U,V,X,Y,XN,YN		
	DOUBLE PRECISION A,TRIANGL
	CHARACTER*17 INPUT
C
C
	OPEN(15,FILE=INPUT,STATUS='OLD')
C
C***	Find starting points of desired areal segments depending on the sector where the
C	outline segment occurs:
C
	IF (ID.EQ.1) THEN				! Case (1): Outline from MA to MB
 5	  READ(15,*) U,V				! Read first point !!!
 	  IF (V.GT.YSTART) THEN				! If first point in file is above 240 pixel line continue searching
	    GOTO 5
	  END IF
C
 	  IF (V.EQ.YSTART) THEN				! If 240 pixel line coincides with first point, set area to zero
	    AREA=0.0					! Initialize AREA for case (1)
	    GOTO 10
	  END IF					
C
	  IF (V.LT.YSTART) THEN				! If first point in file lies below 240 pixel line, 
	    AREA=TRIANGL(XSTART,YSTART,U,V,XF,YF)	! then determine area of first areal increment
	    GOTO 10
	  END IF
	END IF
C
	IF (ID.EQ.2) THEN				! Case (2): Outline segment from MB to MC
 6	  READ(15,*) U,V				! Read first point
 	  IF ((U.EQ.XSTART).AND.(V.EQ.YSTART)) THEN
	    AREA=0.0
	    GOTO 10
	  END IF
	  IF (U.GT.XSTART) THEN
	    AREA=TRIANGL(XSTART,YSTART,U,V,XF,YF)	! Initialize AREA for case (2)
	    GOTO 10
	  END IF
	  GOTO 6
	END IF
C
	IF (ID.EQ.3) THEN
 7	  READ(15,*) U,V
 	  IF ((U.EQ.XSTART).AND.(V.EQ.YSTART)) THEN
	    AREA=0.0
	    GOTO 10
	  END IF
	  IF ((U.EQ.XSTART).AND.(V.GT.YSTART)) THEN
	    AREA=TRIANGL(XSTART,YSTART,U,V,XF,YF)	! Initialize AREA for case (3)
	    GOTO 10
	  END IF
	  GOTO 7
	END IF
C
 	IF (ID.EQ.4) THEN
 8	  READ(15,*,END=999) U,V
 	  IF ((U.EQ.XSTART).AND.(V.EQ.YSTART)) THEN
	    AREA=0.0					! Initialize AREA for case (4)
	    GOTO 10
	  END IF
	  IF ((U.LT.XSTART).AND.(V.EQ.YSTART)) THEN
	    AREA=TRIANGL(XSTART,YSTART,U,V,XF,YF)	! Initialize AREA for case (4)
	    GOTO 10
	  END IF
	  GOTO 8
	END IF
C
 10 	X=U
 	Y=V
 15	READ(15,*,END=999) XN,YN
C
C***	Define stop criteria:
C
	IF (ID.EQ.1) THEN
	  IF (XN.EQ.XEND) THEN
	    GOTO 20
	  ELSE IF (XN.GT.XEND) THEN
	    GOTO 25
	  END IF
C
	ELSE IF (ID.EQ.2) THEN
	  IF (YN.EQ.YEND) THEN
	    GOTO 20
	  ELSE IF (YN.GT.YEND) THEN
	    GOTO 25
	  END IF
C
	ELSE IF (ID.EQ.3) THEN
	  IF (XN.EQ.XEND) THEN
	    GOTO 20
	  ELSE IF (XN.LT.XEND) THEN
	    GOTO 25
	  END IF
C***	  (Note: last case (ID=4) is treated under 999).
	END IF
C
C***	Calculate areal increment and sum up:
C
	A=TRIANGL(X,Y,XN,YN,XF,YF)
	AREA=AREA+A
C
	X=XN					! Permutations for next areal increment
	Y=YN
C
	GOTO 15
C
C***	Make addition of last areal increment for cases (1) through (3):
C
 20	A=TRIANGL(X,Y,XN,YN,XF,YF)
 	AREA=AREA+A
	CLOSE(15)
	RETURN
C
 25	A=TRIANGL(X,Y,XEND,YEND,XF,YF)
 	AREA=AREA+A
 	CLOSE(15)
	RETURN
C
C***	Make addition of last areal increment for case (4)
C	Note: XN and YN will never reach MA because MA is at the beginning of the input file !
C
 999	A=TRIANGL(XEND,YEND,XN,YN,XF,YF)
 	AREA=AREA+A
	CLOSE(15)
	RETURN
C
	END

	    



	
	
	
	

	FUNCTION TRIANGL(U1,V1,U2,V2,U3,V3)
C
C	By Michael Knappertsbusch, 10.2.2003
C
C	Calculates the area of a triangle.
C	Given the cartesian coordinates of three points P1, P2 and P3 e.g. (U1,V1),
C	(U2,V2), and (U3,V3), respectively, this function calulates the area of
C	the triangle defined by P1,P2, and P3.
C	The area TRIANGL is returned to the calling unit in units of the cartesian coordinates.
C	The calculation uses the vector product of Vectors P1P2 and P1P3.
C
	DOUBLE PRECISION U1,V1,U2,V2,U3,V3
	DOUBLE PRECISION TRIANGL
C
	TRIANGL=0.5*ABS((U2-U1)*(V3-V1)-(V2-V1)*(U3-U1))
C
	RETURN
	END










C*******************************************************************************	
	SUBROUTINE EXTREMA(FNAME,UMIN,VA,
     1 UB,VMIN,UMAX,VC,UD,VMAX,I)
C*******************************************************************************
C
C	Modified & corrected 17.2.2003
C	Modified & corrected 28.2.2003
C	Annotated 23.11.2012
C
C	Determines the minima and maxima of U,V-pairs of values and their
C	corresponding U- or V-coordinates in file FNAME.
C
C	Note, that EXTREMA does only determine extreme left and right, upper and lower limits;
C	that are first met in counterclockwise direction. In case of point A values do not
C	necessarily represent the precise apex position or turning positions in this region.
C	In case of B, C and D its always the first extremal point in counterclockwise direction,
C	that is detected and returned to the calling unit.
C	UMIN,UMAX,VMIN,VMAX are returned to the calling unit.
C	EXTREMA also determines the number of points I of the input outline.
C
C	Coordinates of extrema in outline (see sketch):	
C	A(UMIN,VA), C(UMAX,VC)				  D
C	B(UB,VMIN), D(UD,VMAX)				A   C
C							  B
C
C
	INTEGER I
	DOUBLE PRECISION UMIN,VMIN
	DOUBLE PRECISION UMAX,VMAX
	DOUBLE PRECISION U,V,U1,V1
	DOUBLE PRECISION VA,UB,VC,UD
	CHARACTER*17 FNAME
C
	OPEN(15,FILE=FNAME,STATUS='OLD')
C
C	Determination of Umin,Umax,Vmin and Vmax:
C
C	UMIN (point A), and counting the number of points (=I) in the input file:
C
	READ(15,*) UMIN,VA			! Depending on the starting coordinate of the outline
C						! point A is at the beginning or at the end of a vertical line (i.e.
C						! line AA1 in subroutine CORRECT). In the ideal case, point A hits the
C						! y=240 pixels line. If no vertical extension exists in that region
C						! point A represents the apex of the shell.
	I=1					 
 5	READ(15,*,END=100) U,V1
 	I=I+1					! counting the number of points
 	IF (U.LT.UMIN) THEN
	  UMIN=U				! Point A: Here is the first point, that is
C						! hit in the outline on the left side having the
C						! minimum X-coordinate.
	  VA=V1					! YA=corresponding y-coordinate
	END IF 
	GOTO 5
C

 100	REWIND(15)
C
C	UMAX (point C):
C
	READ(15,*) UMAX,VC
 10	READ(15,*,END=200) U,V1
	IF (U.GT.UMAX) THEN
	  UMAX=U				! Point C: Maximum U-coordinate at right side, that is first hit
C						! (coordinates are arranged in counterclockwise direction !)
	  VC=V1					! VC=corresponding y-coordinate
	END IF
	GOTO 10
C
 200	REWIND(15)
C
C	VMIN (point B):
C
	READ(15,*) UB,VMIN
 15	READ(15,*,END=300) U1,V
	IF (V.LT.VMIN) THEN
	  VMIN=V				! Point B: Minimum V-coordinate at the lower end, that is first hit
	  UB=U1					! UB=corresponding U-coordinate
	END IF
	GOTO 15
C
 300	REWIND(15)
C
C	VMAX (point D):
C
	READ(15,*) UD,VMAX
 20	READ(15,*,END=400) U1,V
 	IF (V.GT.VMAX) THEN
	  VMAX=V				! Point D: Maximum V-coordinate at the top, that is first hit
	  UD=U1					! UD=corresponding U-coordinate
	END IF
	GOTO 20
C
 400	CLOSE(15)
C
 	RETURN
 	END







C*******************************************************************************	
	SUBROUTINE KWIDTH10(Y10,XL10,XR10,D10,INPUT)
C*******************************************************************************
C
C	Determines the keelwidths D10 at y-coordinate Y10.
C	Input values are Y10 and INPUT. The variables XL10,XR10 and D10 are returned
C	to the calling unit.
C
	DOUBLE PRECISION X,Y,XN,YN
	DOUBLE PRECISION Y10,D10
	DOUBLE PRECISION XL10,XR10
	CHARACTER*17 INPUT
C
	OPEN(15,FILE=INPUT,STATUS='OLD')
	READ(15,*) X,Y				! Initialize calculations.
 5	READ(15,*,END=100) XN,YN		! Reading next pair of coordinates.
C
C****	Find 10% height x,y coordinates on the left side of the outline profile.
C	Note, that the tracing of the outline occurred in counterclockwise direction
C	and started on the 240 pixel height (e.g. on the top of the spiral side,
C	which was always set at the left hand side).
C
	IF ((Y.GT.Y10).AND.(YN.LT.Y10)) THEN	! Finding x-position of LEFT Y10 height
	  XL10=(X+XN)/2.
	END IF
C
	IF ((Y.LT.Y10).AND.(YN.GT.Y10)) THEN	! Finding x-position of RIGHT Y10 height
	  XR10=(X+XN)/2.0
	END IF
C
C****	Replace current pair of X,Y by XN,YN before reading in next XN,YN values from file:
C
	X=XN
	Y=YN
	GOTO 5
C
 100	D10=XR10-XL10
C
	CLOSE(15)
C
	RETURN
 	END
	
	
	
	
	
	
	
	
C*******************************************************************************	
	SUBROUTINE KWIDTH90(Y90,XL90,XR90,D90,INPUT)
C*******************************************************************************
C
C	Determines the keelwidths D90 at y-coordinate Y90.
C	Input values are Y90 and INPUT. The variables XL90,XR90 and D90 are returned
C	to the calling unit.
C
	DOUBLE PRECISION X,Y,XN,YN
	DOUBLE PRECISION Y90,D90
	DOUBLE PRECISION XL90,XR90
	CHARACTER*17 INPUT
C
	OPEN(15,FILE=INPUT,STATUS='OLD')
	READ(15,*) X,Y				! Initialize calculations.
 5	READ(15,*,END=100) XN,YN		! Reading next pair of coordinates.
C
C****	Find 90% height x,y coordinates on the left side of the outline profile.
C	Note, that the tracing of the outline occurred in counterclockwise direction
C	and started on the 240 pixel height (e.g. on the top of the spiral side,
C	which was always set at the left hand side).
C
C
	IF ((Y.LT.Y90).AND.(YN.GT.Y90)) THEN	! Finding x-position of RIGHT Y90 height
	  XR90=(X+XN)/2.0
	END IF
C
	IF ((Y.GT.Y90).AND.(YN.LT.Y90)) THEN	! Finding x-position of LEFT Y10 height
	  XL90=(X+XN)/2.
	END IF
C
C****	Replace current pair of X,Y by XN,YN before reading in next XN,YN values from file:
C
	X=XN
	Y=YN
	GOTO 5
C
 100	D90=XR90-XL90
C
	CLOSE(15)
C
	RETURN
 	END












C*******************************************************************************
	SUBROUTINE ANGLE(XU,YU,XV,YV,XW,YW,PHI)
C*******************************************************************************
C
C	Calculation of the angle PHI between two vectors VU and VW.
C
C	This subroutine returns positive values of PHI from 0 through 180 degrees.
C	If PHI is greater than 180 degrees it returns the difference
C	of PHI to 360 degrees.
C	Example: If PHI is 270 degrees it returns 90 degrees.
C
C	Input are coordinates of points U,V,W (see sketch below).
C	It is irrelevant, whether vector UV or vector VW is entered first
C	into the subroutine.
C	PHI is returned to the calling unit in degrees.
C
C
C                                              * W(XW,YW)
C		       U(XU,YU) *           .
C                                .        .
C                                 . PHI .
C                                  .  .
C                                   .
C                                V(XV,YV)
C
C
C	The coordinates of U,V and W are:
C	U(XU,YU)	V(XV,YV)	W(XW,YW)
C
C	This subroutine is used to calculate the upper and lower "Keel Angle" of a
C	globorotalid foraminifer in keel view, from which the maximum and
C	minimum extensions of the outline are given.
C	
C
	DOUBLE PRECISION XU,YU,XV,YV,XW,YW
	DOUBLE PRECISION PHI
	DOUBLE PRECISION SUM1,SUM2,SUM
	DOUBLE PRECISION ROOT1,ROOT2
	DOUBLE PRECISION PI,COSPHI
C
	PARAMETER(PI=3.14159)
C
C***	Calculation of PHI:
C
	SUM1=(XU-XV)*(XW-XV)
	SUM2=(YU-YV)*(YW-YV)
	SUM=SUM1+SUM2
C
	ROOT1=SQRT((XU-XV)**2+(YU-YV)**2)
	ROOT2=SQRT((XW-XV)**2+(YW-YV)**2)
C
	COSPHI=SUM/(ROOT1*ROOT2)
	PHI=(180.0/PI)*ACOS(COSPHI)
C
	RETURN
	END












C*******************************************************************************
	SUBROUTINE CORRECT(INPUT,XMIN,YA,XB,
     1 YMIN,XMAX,YC,XD,YMAX,XMA,YMA,XMB,YMB,
     2 XMC,YMC,XMD,YMD)
C*******************************************************************************
C
C	This subroutine is necessary to find the correct keel angle in the
C	axial (keel) view of a globorotalid outline.
C
C	By Michael Knappertsbusch, 11.12.2002.
C	Modified for KeelWidth13 and higher, 23.11.2012.
C
C	The subroutine determines the midpoints of the vertical and horizontal
C	segments in the extreme regions (i.e. left and right, top and bottom,
C	see figure below).
C
C
C
C	   D1(XD1,YMAX)    MD	  D(XD,YMAX)
C                    *.....*.......*
C                   .               .
C                  .                 . outline of a keeled foram in keel position
C   A2(XMIN,YA2)  *                   .
C   A(XMIN,YA) -->* (=240 PIXEL LINE)  .
C                 .                     .
C           MA    *                     * C1(XMAX,YC1)
C                 .                     .
C                 .                     .
C    A1(XMIN,YA1) *                     .
C                  .                    * MC
C                   .                   .
C                    .                  .
C                     .                 .
C                      .                * C(XMAX,YC)
C                       .              .
C                        :             .
C                         .           .
C                          .         .
C               B(XB,YMIN)  *...*...* B1(XB1,YMIN)
C                               MB
C
C	Input to the subroutine are the filename INPUT, and the x,y coordinates of the
C	extreme points A,B,C and D (previously determined by subroutine EXTREMA in main unit).
C	The coordinates of the midpoints MA,MB,MC and MD are returned to the calling unit.
C	(these coordinates are XMA,YMA,XMB,YMB,XMC,YMC, AND XMD,YMD, respectively). From
C	these coordinates the correct keel angles are then calculated in the main unit
C	by using subroutine ANGLE.
C
C	Modifications (added 23.11.2012)
C	In the previous versions up to Keelwith 12.1 the point A2 was not explicitely determined, but point MA
C	was set equal to point A. This was done because outline extraction in program Trace starts at 240 pixel line.
C	In the new version KeelWidth 13.0 and higher versions, the point A2 is explicitely calculated. This affects
C	subsequent determinations of PH1, PHI2 and PHI3.
C
C
C	Requires subroutine MIDPOINT
C
C
C
	DOUBLE PRECISION XMIN,YA,XB,YMIN
	DOUBLE PRECISION XMAX,YC,XD,YMAX
	DOUBLE PRECISION XMA,YMA,XMB,YMB
	DOUBLE PRECISION XMC,YMC,XMD,YMD
	DOUBLE PRECISION YA1,XB1,YC1,XD1
	DOUBLE PRECISION YA2
	DOUBLE PRECISION X,Y
	CHARACTER*17 INPUT
C
C	Open file containing the x,y coordinates of the specimen with the
C	filename INPUT:
C
	OPEN(15,FILE=INPUT,STATUS='OLD')
C
C	Determination of point A1:
C
	YA1=YA					! In case that AA1 is reduced to a singular point
 5	READ(15,*,END=10) X,Y
 	IF ((X.EQ.XMIN).AND.(Y.LT.YA)) THEN
	  YA1=Y					! Point A1 found
	END IF
	GOTO 5
C
 10	REWIND(15)
C
C
C	Determination of point B1
C
 	XB1=XB					! In case that BB1 is reduced to a singular point
 15	READ(15,*,END=20) X,Y
 	IF ((Y.EQ.YMIN).AND.(X.GT.XB)) THEN
	  XB1=X					! Point B1 found
	END IF
	GOTO 15
C
 20	REWIND(15)
C
C
C	Determination of point C1:
C
	YC1=YC					! In case that CC1 is reduced to a singular point
 25	READ(15,*,END=30) X,Y
	IF ((X.EQ.XMAX).AND.(Y.GT.YC)) THEN
	  YC1=Y					! Point C1 found
	END IF
	GOTO 25
C
 30	REWIND(15)
C
C
C	Determination of point D1:
C
	XD1=XD					! In case that DD1 is reduced to a singular point
 35	READ(15,*,END=40) X,Y
 	IF ((Y.EQ.YMAX).AND.(X.LT.XD)) THEN
	  XD1=X					! Point D1 found
	END IF
	GOTO 35
C
 40	REWIND(15)
C
C
C	Determination of point A2 (newly added for KeelWidth13, 23.11.2012):
C
	YA2=YA					! In case that DD1 is reduced to a singular point
 45	READ(15,*,END=50) X,Y
 	IF ((X.EQ.XMIN).AND.(Y.GT.YA)) THEN
	  YA2=Y					! Point A2 found
	  GOTO 50				! Stop searching
	END IF
	GOTO 45
C
 50	CLOSE(15)
C
C


C
C-----------------------------------------------CHECK: CORRECT returns correct coordinates of points A1,B,B1,
C                                                      C,C1,D,D1 and A2. Also the determination of
C						       MA,MB,MC and MD is correct (28.11.2012).
C
C
C
C***	Calculation of midpoint coordinates of distances AA1,BB1,CC1, and DD1 (modified for KeelWidth13, 23.11.2012):
C
	CALL MIDPOINT(XMIN,YA2,XMIN,YA1,XMA,YMA)	! Midpoint MA 
	CALL MIDPOINT(XB,YMIN,XB1,YMIN,XMB,YMB)		! Midpoint MB
	CALL MIDPOINT(XMAX,YC,XMAX,YC1,XMC,YMC)		! Midpoint MC
	CALL MIDPOINT(XD,YMAX,XD1,YMAX,XMD,YMD)		! Midpoint MD
C
	RETURN
	END
	
	
	
	





C*******************************************************************************	
	SUBROUTINE CONVEX(UMA,UMB,UMC,UMD,SPIR,UMBI)
C*******************************************************************************
C
C	Added 23.11.2012, by Michael Knappertsbusch
C
C	Subroutine to calculate the spiral convexity SPIR and umbilical convexity UMBI
C	from the midpoints MA, MB, MC and MD.
C	SPIR is the horizontal distance between the apex of the shell and the vertical midline.
C	UMBI ist the horizontal distance between midline and the rightmost point on the
C	umbilical side. The midline is in the middle between the two parallel vertical lines, that
C	pass through point MD in the upper keel region and through point MB in the lower keel region.
C
C	Note, that only in special cases the uppermost and lowermost points of each keel region
C	consists of a single point. Often, the uppermost and lowermost keelregions consist of
C	several horizontal pixels, that define line D1D in the upper keel region, and line BB1
C	in the lower keel region (see drawing below).
C
C	Input are the coordinates of midpoints MA,MB,MC and MD. The values SPIR and UMBI
C	are returned to the calling unit.
C
C
C                             x
C                          x  x x
C                          x  x x
C                          x  x x
C                          x  x x
C	   D1(UD1,VMAX)    MD(UMD,VMD) 
C                    *.....*..x....* D(UD,VMAX)
C                   .      x  x x   .
C                  .       x  x x    .  outline of a keeled foram in keel position
C   A2(UMIN,VA2)  *        x  x x     .
C   A(UMIN,VA) -->* (=240 PIXEL LINE)  .
C                 .        x  x x       .
C  MA(UMA,VMA)    *        x  x x       * C1(UMAX,VC1)
C                 .___SPIR____x x       .
C                 .        x  x x       .
C    A1(UMIN,VA1) *        x  x x       .
C                  .       x  x x       * MC(UMC,VMC)
C                   .      x  x___UMBI__.
C                    .     x  x x       .
C                     .    x  x x       .
C                      .   x  x x       * C(UMAX,VC)
C                       .  x  x x       .
C                        : x  x x      .
C                         .x  x x     .
C                          .  x x    .
C               B(UB,VMIN) x*.x.*...* B1(UB1,VMIN)
C                          x  x MB(UMB,VMB) 
C                          x  x x       
C                          x  x x       
C                             x
C                             x
C                          midline
C
C
	DOUBLE PRECISION UMA,UMB,UMC,UMD,SPIR,UMBI
	DOUBLE PRECISION UM
C
	UM=(UMB+UMD)/2.0
	SPIR=UM-UMA
	UMBI=UMC-UM
C
	RETURN
	END











	
	
	
	
C*******************************************************************************	
	SUBROUTINE MIDPOINT(XU,YU,XV,YV,XM,YM)
C*******************************************************************************
C
C	Subroutine to calculate the x,y coordinates of the midpoint M of
C	a line extending from point U to V.
C	Coordinates: U(XU,YU), V(XV,YV), M(XM,YM)
C
C	The coordinates XU,YU,XV, and YV enter the subroutine, and
C	XM,YM are returned to the calling unit.
C
	DOUBLE PRECISION XU,YU,XV,YV,XM,YM
C
	XM=(XU+XV)/2.
	YM=(YU+YV)/2.
C
	RETURN
	END
	






C*******************************************************************************
	SUBROUTINE FPOINTS(UMA,VMA,UMB,VMB,
     1  UMC,VMC,UMD,VMD,UF1,VF1,UF2,VF2)
C*******************************************************************************
C
C	By michael Knappertsbusch, 27.1.2003
C
C	This subroutine is necessary for the correct determination of partial areas
C	in subroutine AREA.
C
C	In the ideal case, when the positioning of the specimen on the computer
C	monitor was accurate, the line MDMB is perfectly vertical (see figure below).
C
C	However, this is not always the case, despite it is routinely attempted to
C	(manually) rotate the specimen in vertical position under the microscope. In
C	practice, therefore, there will be a slight deviation in the x-coordinates
C	between points MB and MD.
C
C	If not corrected, this deviation affects the coordinates of the "F-points"
C	F1 and F2 and so calculation of the partial areas in fields (1) through (4)
C	in subroutine AREA might become errorneous.
C
C	In subroutine FPOINTS the coordinates of the "F-points" F1 and F2 are calculated
C	precisely, so that the error in the area computation in subroutine AREA vanishes.
C
C
C
C
C	   D1(XD1,YMAX)        MD  D(XD,YMAX)
C                          *....*...*
C                    .....      |   .
C                  ..     (4)   |    .   Outline of a keeled foram in keel position
C   A2(XMIN,YA2)  *             |     .
C   A(XMIN,YA) -->* (=240 PIXEL LINE)  .
C                 .             |       .
C           MA    *-------------+ F1    * C1(XMAX,YC1)
C                 .             |       .
C                 .             |  (3)  .
C    A1(XMIN,YA1) *    (1)      |       .
C                  .         F2 +-------* MC
C                   .           |       .
C                    .          |       .
C                     .         |       .
C                      .        |  (2)  * C(XMAX,YC)
C                       .       |      .
C                        :      |      .
C                         .     |     .
C                          .    |    .
C               B(XB,YMIN)  *...*...* B1(XB1,YMIN)
C                               MB
C
C
	DOUBLE PRECISION UMA,VMA,UMB,VMB
	DOUBLE PRECISION UMC,VMC,UMD,VMD			! Coord's of MA,MB,MC and MD.
	DOUBLE PRECISION UF1,VF1,UF2,VF2			! Coord's of "F-points" F1 and F2.
C
	IF (UMB.EQ.UMD) THEN					! Line MBMD is perfectly vertical.
	  UF1=UMB
	  VF1=VMA
C
	  UF2=UMB
	  VF2=VMC
	ELSE
	  UF1=(VMA-VMB)*(UMD-UMB)/
     1        (VMD-VMB)+UMB
	  VF1=VMA
C
	  UF2=(VMC-VMB)*(UMD-UMB)/
     2        (VMD-VMB)+UMB
	  VF2=VMC
C
	END IF
C
	RETURN
	END









