6. Appendix - MorphCol program listings

6.1. Listing of program Trace33_batch.f

PROGRAM TRACE
C
C !!! Processing in batch mode. !!!
C Version 1.2.3, April 28, 1997 by M. Knappertsbusch,
C modified to batch mode on 21.12.1999
C modified to write errorneous images to external file. 21.1.2000
C modified to record missing images, that are read from FILE_LIST. 24.1.2000
C modified by changing MAG from INTEGER to REAL. 18.12.2000
C modified by changing the calibration (New XPREC and YPREC). 22.12.2000
C
C Program to trace the outline of a series of digitized light images
C in batch mode. Input images are grey-level pictures with 480 lines
C and 680 pixels per line.
C
C
C !!!!! Output are the X,Y coordinates of the outline
C in COUNTER-CLOCKWISE direction !!!!!
C Program includes conversion of pixels in micrometers.
C Program calibrated to the Leica MZ6 binocular stereo microscope.
C
C Input files are:
C - A File with name FILE_LIST containing a list of all
C names of the input images (filename) to be traced,
C and the magnifications, that were applied to each image.
C FILE_LIST is a character variable and can bear any name
C up to 20 characters in length.
C
C - All image files to be traced. The names of these raw images
C are interpreted through the character variable INPUT.
C INPUT is exactly 5 characters long.
C
C - The program calls for the sample name through the
C character variable SAMPLE. SAMPLE is esactly 11 characters
C long (including orientation codes K or U for Keel or Umbilical).
C
C INPUT and SAMPLE are used to compose the name of the corresponding
C output file (suffix _T).
C
C Output files are:
C - One file OUTPUT for every traced sample, containing the x,y
C coordinates of the outline. The character variable OUTPUT is
C 17 characters long (including the suffix _T).
C
C Note: All file handling for morphometric analyses follows a name convention:
C
C filename = Name of sample, 15 characters long.
C filename_T = Input file for raw curves with cartesian coordinates, 17 characters long.
C filename_POL = Polar coordinates of filename, 19 characters long.
C filename_INT = Interpolated (cartesian) coordinates of filename, 19 characters long.
C filename_NORM = Normalized curve, 20 characters long.
C
C
C*** Note: MAG must be real, otherwise results may be wrong !
C
INTEGER I,N,N0,X,Y,THRES,NPTS,NF
INTEGER GREY(1:8), K(1:8), A(1:8,1:8)
LOGICAL IMAGEX
REAL MAG
DOUBLE PRECISION XPREC,YPREC,XRES,YRES
CHARACTER*1 CHAR(1:8),ANSWER,COMMA
CHARACTER*5 INPUT
CHARACTER*17 OUTPUT
CHARACTER*11 SAMPLE
CHARACTER*20 FILE_LIST
PARAMETER (THRES=49)
C
C In A(I,J) is the local numbering of points:
C

DATA ((A(I,J),I=1,8),J=1,8)/1,
* 2,3,4,5,6,7,8,
* 2,3,4,5,6,7,8,1,
* 3,4,5,6,7,8,1,2,
* 4,5,6,7,8,1,2,3,
* 5,6,7,8,1,2,3,4,
* 6,7,8,1,2,3,4,5,
* 7,8,1,2,3,4,5,6,
* 8,1,2,3,4,5,6,7/
C
C Experimental setup and calibration:
C
C A Leica MZ6 Stereo-Binocular eqipped with a CF 11/2 colour
C CCD-camera 795(H)x596(V) pixels from Kappa, and imaged using Nih-Image 1.61
C on a Power Macintosh 8500 computer. The camera is connected to the
C microscope via a 1x Cmount.
C The magnification to be entered is calculated as the product
C of the objective lens inserted (1x or 2x) times the number indicated
C at the zoom-knob. Example: if the objective lens is 1x and the
C number at the zoom-knob indicates 4.0 then the magnification to be
C entered is 4.0.
C
C In the following factors are calculated to convert pixels into µm with the
C CCD camera, using the NTSC (640x480 pixels format in NIH):
C (XPREC and XRES are in horizontal, YPREC and YRES are in vertical direction).
C This calibration includes now a newer one, which was done with the 0.01mm scale
C (=small scale) and the 0.1mm scale (large scale from Leica). This has not much
C changes in the final result but the linearity of XPREC and YPREC over the entire
C range of magnifications from 0.63x to 8x could so be demonstrated (see tests
C from 21.12.2000).
C
C XPREC = 0.11984 * MAG + 0.00048188 (in pixel/µm)
C YPREC = 0.11926 * MAG + 0.00074223 (in pixel/µm)
C
C The necessary resolutions XRES and YRES are the calculated as
C
C XRES = 1/XPREC (in µm/pixel)
C YRES = 1/YPREC (in µm/pixel)
C
C Specifications on computer side:
C Frame grabber: built in Quicktime framegrabber.
C Video control: Grey-scale, full-size, NTSC format
C Digital image: resolution: 72 pixel per inch, image size: 640x480 pixel
C = Quicktime format).
C Pixel aspect ratio on computer screen: 1
C
C
WRITE(9,*) ' Program Trace'
WRITE(9,*) 'Version 3.3 calibration from 21.12.2000'
WRITE(9,*) ' Processing in batch mode'
WRITE(9,*) ' '
WRITE(9,*) '. . .Enter sample name (11 chars, with K or U). . .'
READ(9,102) SAMPLE
102 FORMAT(A11)
C
WRITE(9,*) '. . .Output in Pixels (1) or Microns (2) ?. . .'
READ(9,4) ANSWER
4 FORMAT(A1)
IF (ANSWER.EQ.'1') THEN
WRITE(9,*) 'Output is in pixels'
ELSE
WRITE(9,*) 'Output is in µm'
WRITE(9,*) ' '
END IF
C
C Get filenames of images and magnifications from FILE_LIST, and start opening
C and reading for first image:
C
WRITE(9,*) '. . .Enter list containing image
1 names (20 chars) and MAG. . .'
WRITE(9,*) 'Format: 5 chars, comma ,1 real value (2 decimals)'
WRITE(9,*) 'Note: MAG=Objective lens x zoom factor'
WRITE(9,*) 'Example: 0101r,3.20 (specimen=0101r,
2 magnification=3.20x)'
READ(9,101) FILE_LIST
101 FORMAT(A20)
C
C Initialization of NF (NF= number of files traced):
C
NF=0
C
C Recording errors during program run:
C
OPEN(21,FILE='Image_Errors',STATUS='NEW')
WRITE(21,96) SAMPLE
96 FORMAT('Errors in Images of Sample: ',A11)
OPEN(17,FILE='Missing_images',STATUS='NEW')
C
C Now reading filenames etc:
C
OPEN(20,FILE=FILE_LIST,STATUS='OLD')
100 READ(20,110,END=999) INPUT,COMMA,MAG
110 FORMAT(A5,A1,F4.2)
C
C Check first, whether image with name INPUT exists or not. If it does not exist:
C print a message and continue reading next image.
C
INQUIRE(FILE=INPUT,EXIST=IMAGEX)
IF (.NOT.IMAGEX) THEN
WRITE(9,*) 'Error: image ',INPUT,' is missing'
WRITE(17,97) INPUT
97 FORMAT(A5)
WRITE(9,*) '. . .Program continues. . .'
GOTO 100
END IF
C
C Counting the number of images (=NF) treated:
C
NF=NF+1
C
C Initialization of NPTS: (=Number of points in outline):
C
NPTS=1
C
C Determination of the names of the OUTPUT files (17 chars long), and
C writing information to screen:
C
OUTPUT=(SAMPLE(1:11)//INPUT(1:4))//'_T'
WRITE(9,111) OUTPUT,MAG
111 FORMAT(' Specimen: ',A17,' at magnification: ',F4.2,'x')
C
C
C Determination of XRES and YRES (has been calibrated using the NTSC
C (640x480 pixel) format in NIH-Image:
C
XPREC=0.11984*MAG+0.00048188
YPREC=0.11926*MAG+0.00074223
C
XRES=1/XPREC
YRES=1/YPREC
C
C
C WRITE(9,*) '. . .calculating. . .'
C
OPEN(15,FILE=INPUT,ACCESS='DIRECT',
* RECL=1,FORM='FORMATTED',STATUS='OLD')
OPEN(16,FILE=OUTPUT,STATUS='NEW')
C
C Find first point of outline in the middle of the
C TV-screen (this is line 240; its record is calculated
C as 152960=(240-1)*640):
C
DO 1, I=1,640
N=152960+I
READ(15,50,REC=N) CHAR(1)
GREY(1)=ICHAR(CHAR(1))
IF (GREY(1).GE.THRES) THEN
N0=N
GOTO 5
END IF
1 CONTINUE
C
5 CALL COORD(N,X,Y)
C WRITE(9,*) X,Y
C
C Output of results to file:
C
IF (ANSWER.EQ.'1') THEN
C Output of X and Y in pixels:
WRITE(16,*) X,',',Y
ELSE
C Output of X and Y in micrometers:
WRITE(16,*) X*XRES,',',Y*YRES
END IF
C
C
C STOP Criterium if there is an error in the image:
C (NPTS indicates the number of points written to the outline file; usually,
C an outline of forams at a magnification of <4x has less than 2500 points).
C The errorneous INPUT images are written into external file with name
C "Image_Errors", that can be used as input for subsequent runs.
C
NPTS=NPTS+1
IF (NPTS.GT.2500) THEN
WRITE(9,119) OUTPUT
WRITE(21,110) INPUT,COMMA,MAG
119 FORMAT('Error in image ',A17,'. Reading next image')
GOTO 100
END IF
C
C
C
C Local initialization: At each point (pixel) of the outline (record number N)
C the eight surrounding pixels are read, then their grey-levels
C are determined, and the pixels are numbered from 1 to 8:
C
C
C 2 3 4
C
C 1 N 5
C
C 8 7 6
C
C
READ(15,50,REC=N-1) CHAR(1)
READ(15,50,REC=N-641) CHAR(2)
READ(15,50,REC=N-640) CHAR(3)
READ(15,50,REC=N-639) CHAR(4)
READ(15,50,REC=N+1) CHAR(5)
READ(15,50,REC=N+641) CHAR(6)
READ(15,50,REC=N+640) CHAR(7)
READ(15,50,REC=N+639) CHAR(8)
50 FORMAT(A1)
C
C Determination of the grey-levels:
C
DO 10, I=1,8
GREY(I)=ICHAR(CHAR(I))
10 CONTINUE
C
C Assignment of the local pixel number (i.e. 1 through 8)
C to the record number N for each pixel.
C (i.e. in K are the record numers of the pixels 1 to 8 stored,
C with K having the idices from i to 8):
C
K(1)=N-1
K(2)=N-641
K(3)=N-640
K(4)=N-639
K(5)=N+1
K(6)=N+641
K(7)=N+640
K(8)=N+639
C
C Search algoritm for the next point of the outline: The algoritm
C searches first for the local point to start with: This first
C point has the grey-level > Threshold. From this initial
C local point it searches the next point, where the greylevel is
C lower than Threshold. The pixel, which is previous to that
C point is the next point of the outline (i.e. the next center for
C the search algoritm).
C
C Look for first point of lokal environment-search:
C
DO 20, I=1,8
IF (GREY(I).GE.THRES) THEN
C
C First point of search path found. Check every next surrounding point in
C clockwise direction:
C
DO 30, J=1,8
IF (GREY(A(I,J)).LT.THRES) THEN
IF (A(I,J).EQ.1) THEN
N=K(8)
ELSE
N=K(A(I,J)-1)
END IF
C
C Stop-Criterium:
C
IF (N.EQ.N0) THEN
C WRITE(9,*) 'Outline closed'
CLOSE(16)
CLOSE(15)
C
C Reading next file:
GOTO 100
C
END IF
C
GOTO 5
END IF
30 CONTINUE
END IF
20 CONTINUE
C
999 WRITE(9,*) ' '
WRITE(9,*) NF,' images traced and written to files'
WRITE(9,*) 'One file "Image_Errors" written'
WRITE(9,*) 'One file "Missing_images" written'
WRITE(9,*) 'Enter a key to end the program'
PAUSE 888
888 CONTINUE
C
STOP
END

SUBROUTINE COORD(M,U,V)
C
C Determines the cartesian coordinates U and V as a function of
C the record number M representing a pixel in the image.
C The image has 480 lines and 640 pixels per line.
C The record number starts with M=1 (representing the first
C pixel in line 1).
C
INTEGER M,U,V
C
C This determination of U and V is valid for all pixels except those
C that are at the end of each line (i.e. pixel no. 640, 1280,...).
C
U=MOD(M,640)
V=479-INT(M/640)
C
RETURN
END

6.2. Listing of program XYPlot2.f

GLOBAL DEFINE
INCLUDE "Windows.inc"
END

PROGRAM XYPlot
C
C 30.10.2003 by Michael Knappertsbusch
C
C Program to draw X,Y plots from external X,Y data to the screen.
C
INTEGER I
REAL X,Y
CHARACTER*20 LIST,FILENAME
PARAMETER(RHO=1) ! Radius of circles
C
RECORD /Rect/rrect ! Define rectangle structure rrect
C
rrect.top=0 ! Size of rrect
rrect.left=0
rrect.bottom=600
rrect.right=600
C
C
WRITE(9,*) '***********************************'
WRITE(9,*) '* *'
WRITE(9,*) '* Program XYPlot, V.2 (Batch mode) *'
WRITE(9,*) '* *'
WRITE(9,*) '* by Michael. Knappertsbusch, 30.10.2003 *'
WRITE(9,*) '* Draws X,Y data to screen at intervals of *'
WRITE(9,*) '* 5 outlines *'
WRITE(9,*) '* *'
WRITE(9,*) '***********************************'
WRITE(9,*) ' '
WRITE(9,*) '. . .Enter filename with list of files'
WRITE(9,*) 'that contain X,Y data. . .'
WRITE(9,*) ' (max 20 chars long)'
READ(9,5) LIST
5 FORMAT(A20)
C
C
C Main loop:
C
CALL EraseRect(rrect) ! Clear screen
OPEN(14,FILE=LIST,STATUS='OLD')
C
I=0 ! Initialize counter I
10 K=0
I=I+1
READ(14,5,END=999) FILENAME ! Reading the filenames from the list
OPEN(15,FILE=FILENAME,STATUS='OLD')
20 K=K+1
IF ((MOD(I,5).EQ.0).AND.(K.EQ.1)) THEN ! Avoid plotting first point if MOD(I,5) is zero
GOTO 30
END IF
READ(15,*,END=100) X,Y ! Reading X,Y coordinates from a file
CALL DrawCircle(RHO,100+X/5,100+Y/5) ! Plot points (=circles of radius 1)
C
C Draw intervals of 5 outlines and then pause:
C
30 IF ((MOD(I,5).EQ.0).AND.(K.EQ.1)) THEN
PAUSE
CALL EraseRect(rrect) ! Erase previous 5 outlines from screen
END IF
C
GOTO 20
C
100 CLOSE(15)
GOTO 10
999 CLOSE(14)
PAUSE
C
STOP
END
SUBROUTINE DrawCircle(Radius,Xm,Ym)
C
C Draws a circle with radius Radius, and Center with coordinates
C at Xm,Ym.
C
RECORD /Rect/ r
RECORD /WindowRecord/ front
POINTER (p_front, front)
C
p_front=FrontWindow()
CALL SetPort(VAL4(p_front))
C
r.top=Ym-Radius ! Defines a rectangle structure
r.left=Xm-Radius
r.bottom=Ym+Radius
r.right=Xm+Radius
C
CALL FrameOval(r) ! Draws the oval (or circle)
C
RETURN
END

6.3. Listing of program Sprep53.f

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
C 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 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=2*PI/3
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

6.4. Listing of program KeelWidth93.f

program KeelWidth
C
C Version 9.3 by Michael Knappertsbusch
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
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. The _T files are read and treated sequentially.
C
C Output files:
C There is one 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 Example:
C reading the specimen 502A106125K0001_T from "MEASUREMENTS" and getting menardA as
C species code, 6chambers in the last whorl and a sinistral coiling direction, the output
C filename, which contains the results will be RESKW_502A106125K_6Ch_s.
C
C
C List of subroutines:
C EXTREMA, KWIDTH10, KWIDTH90, ANGLE, CORRECT, MIDPOINT, FPOINTS, ARSEG
C
C List of functions:
C TRIANGL
C
C
INTEGER N,NPOINTS,LAN
LOGICAL FILEX
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 ! Q is 0.1 for D10 at 10% and D90 at 90% of the vertical extension.
DOUBLE PRECISION XMIN,YMIN,XMAX,YMAX ! XMIN,XMAX are the leftmost or rightmost X-
C ! coordinates of the outline,
C ! YMIN,YMAX are the lowermost or uppermost Y-coordinates of the
C ! outline.
DOUBLE PRECISION YA,XB,YC,XD ! X, and Y-coordinates of the corresponding values for the
C ! extrema (see in subroutine EXTREMA).
DOUBLE PRECISION Y10,Y90,D10,D90 ! Y10,Y90 are the y-coordinates at the lower or upper
C ! keelwidths, D10
DOUBLE PRECISION D9010,DMIN,DMAX ! and D90 are the values of the lower and upper keelwidths,
C ! D9010 is
DOUBLE PRECISION DX,DY ! the absolute difference between D90 and D10, DMIN is the
C ! smaller and DMAX
DOUBLE PRECISION XL10,XR10,XL90,XR90 ! DMAX is the larger of D10 and D90, and DX is the difference
C ! 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
C ! 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
C ! areas".
DOUBLE PRECISION AREA3,AREA4
DOUBLE PRECISION XSI1,XSI2,XSI3,XSI4 ! Areas of triangles (1) to (4) for determination of "keel
C ! 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:
C ! Midpoint between B and B1;
DOUBLE PRECISION XMC,YMC,XMD,YMD ! XMC,YMC: Midpoint coords between C and C1; XMD,YMD:
C ! 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
C degrees.
DOUBLE PRECISION TRIANGL ! Function TRIANGL, that determines the area of the triangle
C ! given by 3 points.
DOUBLE PRECISION DEV,SUM1 ! DEV: Control variable (in %) for estimation of deviation of
C ! sum of (AREA1+AREA2+AREA3+AREA4)
C ! from AREA (given in file "MEASUREMENTS" during determination of
C ! ARESLO and ARESUP.
C
CHARACTER*66 HEADER ! Header line in file MEASUREMENTS
CHARACTER*7 SPECIES ! Species name (same name encoding as used in programs
C ! KeelBend, Homolo, etc).
CHARACTER*17 INPUT ! Name of the file with the outline coordinates (_T files, 17
C ! chars long).
CHARACTER*31 OUTPUT ! Name of the file with the results.
CHARACTER*1 COIL ! Coiling direction.
CHARACTER*2 NCHAMB ! NCHAMB: Number of chambers.
C
PARAMETER (Q=0.1) ! Modify here to change the Y-positions of D10 and D90.
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 9.3 By Michael Knappertsbusch *'
WRITE(9,*) '* *'
WRITE(9,*) '**************************************'
WRITE(9,*) ' '
WRITE(9,*) ' '
C
C Recording errors during program run:
C
OPEN(25,FILE='Errors_during_KeelWidth',STATUS='NEW')
C
WRITE(9,*) 'Notes:'
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
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), separated by commas'
WRITE(9,*) 'Example: menardA, 6,s'
READ(9,14) SPECIES,NCHAMB,COIL
14 FORMAT(A7,1X,A2,1X,A1)
C
WRITE(9,*) ' '
WRITE(9,*) ' '
WRITE(9,2)
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 %)')
C
C
C Now, compose the name of the output file:
C
LAN=2 ! Sets the value for NCHAMB to 2 characters
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
OPEN(16,FILE=OUTPUT,STATUS='NEW') ! File RESKW... contains the results
WRITE(16,*) 'Notes:' ! Write various notes to RESKW...
WRITE(16,15) Q*100.0,100*(1.0-Q)
WRITE(16,1)
WRITE(16,*) ' '
WRITE(16,*) ' '
WRITE(16,2) ! Labelling the columns in RESKW...
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
DX=XMAX-XMIN ! DX is the maximum extension of the shell in X direction
C ! (=spiral height X).
DY=YMAX-YMIN ! DY is the maximum extension of the shell in Y direction
C ! (=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 (default: Q=0.1)
Y90=YMIN+(1.0-Q)*(YMAX-YMIN) ! Y-position of upper keelwidth
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
C ! profile)
CALL KWIDTH90(Y90,XL90,XR90,D90,INPUT) ! Upper keelwidth (at RATIO (default 90)% of height of
C ! 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
CALL CORRECT(INPUT,XMIN,YA,XB,
1 YMIN,XMAX,YC,XD,YMAX,XMA,YMA,XMB,YMB,
2 XMC,YMC,XMD,YMD)
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 PHI1 and PHI2 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
C
C
C
C**** Calculation of "keel residual area" (added 30.1.2003):
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
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

WRITE(16,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
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)
C
C
C**** Reading next file:
C
GOTO 5
C
999 WRITE(9,*) ' '
WRITE(9,*) 'All files treated'
PAUSE1000
1000 CONTINUE
C
STOP
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
C ! continue searching
GOTO 5
END IF
C
IF (V.EQ.YSTART) THEN ! If 240 pixel line coincides with first point,
C ! 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
C Determines the minima and maxima of U,V-pairs of values and their
C corresponding U- or V-coordinates in file FNAME.
C UMIN,UMAX,VMIN,VMAX are returned to the calling unit.
C It 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 ! In the ideal case point A coincides with the
I=1 ! first data point in file INPUT (e.g. line Y=240 pixels)
5 READ(15,*,END=100) U,V1
I=I+1 ! counting the number of points
IF (U.LT.UMIN) THEN
UMIN=U ! Point A: Minimum x-coordinate (at left side)
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)
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)
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)
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 The coordinates of the points U,V,W enter the subroutine, and
C the angle PHI is returned to the calling unit in degrees (see
C sketch below).
C U W
C V
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 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
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 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 Important note: Point A2 is not determined in this subroutine; point MA is set
C equal to point A, because outline extraction in program Trace starts at 240 pixel line.
C
C
C Requires subroutine MIDPOINT
C
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 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
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
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
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
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 CLOSE(15)
C
C
C
C*** Calculation of midpoint coordinates of distances AA1,BB1,CC1, and DD1:
C
XMA=XMIN ! Choose point MA equal to point A, because outline
YMA=YA ! extraction in program Trace starts at 240 pixel
C ! line.
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 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 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

6.5. Listing of Progra, KeelBend7.f

PROGRAM KEELBEND
C Version 7
C
C Program to determine bending of outline in upper and lower keel region.
C
C By Michael Knappertsbusch, 22 August 2000.
C
C Modified 29 August 2000 (addition of optional file output of x,y coordinates
C of keel region)
C Modified 5 September 2000 (Age in Results).
C Modified 3 January 2001 (reading sample names from file MEASUREMENTS,
C write them to file RESULTS and also write the bending characteristicts
C and center coordinates of Kruemmungskreise to that file.
C Modified 5 January 2001: Addition of optional regression of order 3 to
C approximate the keel region by a cubic spline, and appropriate modification
C of subroutine BEND and SLE (variable dimensioning in SLE).
C Version 6 combines file handling for entry of _T, _INT, or _NORM files in
C the same program.
C Modified 10 May 2001: Addition of optional output for extended result file. The program
C asks whether the results shall be written to file RESULT (option 1, simple output as
C in version 6) or, whether the program shall write the results to a file
C RES_samplename_species_#Chambers_coiling (option 2, extended output). The name of
C this extended fil is calculated interactively and from the respective name of the
C sample. In option 2 additional variables are written such as R=X/Y, species names, the
C # of chambers in the last whorl, and the coiling direction. This file can then further be used
C for Program Compose.f (preparation of plots of individual measurements), or for
C Program Stats4.f (calculation of statistics and plot preparation in batch mode).
C
C
C !!! Batch processing mode !!!
C
C
C INPUT:
C Given is the closed plane curve represented by X,Y coordinates
C (i.e. an outline of a microfossil). That Data points are
C stored in a file with name of character value of INPUT.
C The names of the INPUT files are listed, together with measured variables,
C from file MEASUREMENTS. MEASUREMENTS was generated by program Sprep52.f
C and contains variables like the Age of the sample (AGE), the number of
C points in the original outline (N), the spiral height of the shell (X, in µm),
C the width or aequatorial diameter of the shell (Y, in µm), and the
C surface area in keel view (in mm2).
C
C This program searches for a suite of points in the vicinity
C of the keel area in keel view of a globorotalid foraminifer. These
C points can then be used to measure the keel development (in keel
C view) by approximation of the keel outline segment by a fifth order
C polynom and determination of bending around the keel.
C
C The keel area is located by recognizing the point with the maximum and
C minimum Y coordinates, the point with the second maximum and minimum Y
C coordinates, the point with the third maximum and minimum
C Y coordinates, etc., until the inferior points with number K.
C The program exports E points (note, that in a closed curve several
C points can attain the same Y value, so that often E > K.
C The case of E=K is only fulfilled if every point has a different Y value.
C
C OUTPUT:
C The program generates four types of files:
C For every input file there are two files written containing the X,Y
C coordinates of the upper or lower keel region, respectively.
C The output of these files is optional.
C Their names follow the name convention for all shape analysis programs
C (see subroutine CURVREG5) and they are stored in the character
C variable OUTPUT in subroutines MAXIMA or MINIMA.
C The variables for each specimen are stored in a third file type, which
C is either called RESULT (=option 1) or RES_sample_Spec_nCh_coil (=option 2),
C depending on the option selected. The actual specimen to be treated is
C read from file MEASUREMENTS.
C In the fourth file type called Missing_files, the program records
C input files, that are listed in MEASUREMENTS, but that were not found
C during program run.
C
C
C
C The maxima and minima and subordinate extrema values of the Y
C coordinates are stored in vectors YMAX and YMIN. The absolute maxima
C and minima are contained in the first vector components, the second
C extrema values are in the second vector component, and so forth.
C The number of (sub)extrema desired is given by E.
C
C The coordinates of individual points of the curve are read in
C as X and Y.
C
C
C File name handling convention:
C All file names for morphometric analyses follow the same name convention:
C
C filename = Name of sample; 15 chars long.
C filename_T = File of raw curve with X,Y coordinates; 17 chars long.
C filename_POL = File with polar coordinates; 19 chars long.
C filename_INT = File with (cartesian) interpolated coords, 19 chars long.
C filename_MaxK = File with coords of upper keel region, 20 chars long.
C filename_MinK = File with coords of lower keel region, 20 chars long.
C filename_NORM = File with normalized curve, 20 chars long.
C RES_sample_species_nCh_coil = Extended output file with variables, 29 chars long;
C sample is 11 chars long, species is 7 chars long,
C nCh is 3 chars long, coil is 1 char long.
C
C
C Subroutines, that are called are:
C MAXIMA
C MINIMA
C CREG3 (Modified from CURVEREG3)
C CREG5 (Modified from CURVEREG5)
C BEND3
C BEND5
C
C
INTEGER E,NMX,NMI,K,N,Q ! N is # of points in original traced outline (_T file). Q is # of
C ! points in file INPUT
INTEGER ANSWER,TYPE,DEGREE ! DEGREE is 3 for cubic spline or 5 for 5th order polynomial regression.
DOUBLE PRECISION X1,YMAX,X2,YMIN,Y1,Y2 !Y1 and Y2 are the y coordinates of the centers
C of the upper and lower Kruemmungskreis, respectively
DOUBLE PRECISION RHOMAX,RHOMIN,KMAX,KMIN
DOUBLE PRECISION DX,DY,R,AREA ! DX=XMAX-XMIN in µm, DY=YMAX-YMIN in µm, AREA=Area
C ! in mm2
REAL AGE ! R=DX/DY (=flatness of shell in keel view)
DIMENSION A(6),B(6)
LOGICAL FILEX
CHARACTER*1 ANS,NCHAMB,COIL ! ANS=simple or extended output,NCHAMB=# Chambers,
C ! COIL=coiling direction.
CHARACTER*7 SPECIES ! Name of the actual species.
CHARACTER*5 ERR ! ERR indicates if the Determinant in SLE is zero.
CHARACTER*20 INPUTMAX,INPUTMIN
CHARACTER*15 SPECIMEN
CHARACTER*25 TINTNORM ! Qualifier for traced-, normalized-, or interpolated files.
CHARACTER*20 INPUT ! Maximum length of ALL filenames is 20 characters
CHARACTER*29 OUTPUT ! Contains the name of the extended output file.
CHARACTER*66 HEADER
C
C
WRITE(9,*) 'Program to find bending of outline in keel region of a
* globorotalid foraminifera,'
WRITE(9,*) 'that is oriented in keel view'
WRITE(9,*) ' '
WRITE(9,*) 'Batch processing mode'
WRITE(9,*) ' '
WRITE(9,*) '. . .Reading MEASUREMENTS. . .'
C
C Open existing file MEASUREMENTS, read sample age header:
C
OPEN(27,FILE='MEASUREMENTS',STATUS='OLD')
READ(27,*) AGE
WRITE(9,*) 'Sample age is ',AGE,'( Ma)'
READ(27,4) HEADER
4 FORMAT(A66)
C
C Prepare for determination of name of the output file:
C
READ(27,2) SPECIMEN
2 FORMAT(A15)
REWIND(27) ! (The following 2 lines necessary because of statement label 8)
READ(27,*) AGE
READ(27,4) HEADER ! (Sets cursor to line 3)
C
C Now specify which type of input data is desired:
C
WRITE(9,*) ' '
WRITE(9,*) '_T, _INT, or _NORM files ?'
WRITE(9,*) ' (enter 1,2 or 3)'
READ(9,*) TYPE
C
C Determine qualifier for input files:
C
IF (TYPE.EQ.1) THEN
TINTNORM='traced (_T) files'
WRITE(9,3) TINTNORM
ELSE IF (TYPE.EQ.2) THEN
TINTNORM='interpolated (_INT) files'
WRITE(9,3) TINTNORM
ELSE IF (TYPE.EQ.3) THEN
TINTNORM='normalized (_NORM) files'
WRITE(9,3) TINTNORM
END IF
3 FORMAT('. . .Analyzing ',A25,'. . .')
C
C Now read info how the keel region is defined and select degree of polynomial regression
C for keel approximation:
C
WRITE(9,*) ' '
WRITE(9,*) 'Now enter the number of points to be'
WRITE(9,*) 'included for defining the keel region'
WRITE(9,*) ' '
WRITE(9,*) 'Enter an integer value (example: 50)'
READ(9,*) E
WRITE(9,*) ' '
WRITE(9,*) 'Indicate degree of polynomial regression
* in keel region approximation (enter 3 or 5)'
WRITE(9,*) '3 = cubic spline, 5 = 5th order
* polynomial regression'
READ(9,*) DEGREE
C
IF (DEGREE.EQ.3) THEN
WRITE(9,*) 'Cubic spline approximation in
* keel region'
ELSE IF (DEGREE.EQ.5) THEN
WRITE(9,*) '5th order polynomial regression in
* keel region'
END IF
C
C
C Decide whether the output file is simple (RESULT, = option 1) or extended for
C subsequent treatment with programs Compose.f and Stats4.f (= option 2):
C
WRITE(9,*) 'Simple output of variables to file RESULT (=1)'
WRITE(9,*) 'or extended output (=2;
1 for Compose.f and Stats4.f) ?'
WRITE(9,*) 'Enter 1 or 2'
READ(9,'(A1)') ANS
C
IF (ANS.EQ.'1') THEN ! Simple output to file RESULT
C
C Open outputfile RESULT.
C
OPEN(18,FILE='RESULT',STATUS='NEW')
WRITE(18,*) 'Measurements and curvature in keel region:'
WRITE(18,*) ' '
WRITE(18,6) TINTNORM
6 FORMAT('Specimen outlines are ',A25)
WRITE(18,7) E,AGE,DEGREE
WRITE(18,*) ' '
WRITE(18,*) 'Notes:'
WRITE(18,*) '#_Tpts indicates the number of points in the
* traced (_T) file'
WRITE(18,*) '#_Acpts indicates the number of points in the
* actual input file'
WRITE(18,*) 'X1,Y1 and X2,Y2 are coordinates of
* centers of upper and lower Kruemmungskreise, respectively,
* of the ',TINTNORM
WRITE(18,*) 'Rupkeel and Rlokeel are the absolute values of
* "upper and lower Kruemmungsradii" (in µm), respectively,
* of the ',TINTNORM
WRITE(18,*) 'Bupkeel and Blowkeel are the upper and lower
* "Kruemmung", respectively of the ',TINTNORM
7 FORMAT('No. of keel coordinates analyzed: ',I4,
*', Age = ',F5.2,' Ma, degree of polynomial regression: ',I2)
WRITE(18,*) ' '
WRITE(18,*) 'Sample, Age (Ma) #_Tpts #_Acpts X (µm)
* Y (µm) Ar(mm2) Rupkeel Rlokeel Bupkeel Blowkeel
* X1 Y1 X2 Y2'
C
ELSE IF (ANS.EQ.'2') THEN ! Extended output...
C
C Determination of name of output file:
C
WRITE(9,*) 'Enter species name (7chars),
* # of chambers (1 char), and'
WRITE(9,*) 'the coiling (s or d), separtated by commas'
WRITE(9,*) 'Example: menardA,6,s'
READ(9,*) SPECIES,NCHAMB,COIL
C
C Now, compose the name of the output file:
C
OUTPUT='RES_'//SPECIMEN(1:11)//'_'
*//SPECIES//'_'//NCHAMB//'Ch_'//COIL
C
OPEN(18,FILE=OUTPUT,STATUS='NEW')
C
WRITE(18,*) 'Extended measurements
1 and curvature in keel region:'
WRITE(18,*) ' '
WRITE(18,6) TINTNORM
WRITE(18,7) E,AGE,DEGREE
WRITE(18,*) ' '
WRITE(18,*) 'Notes:'
WRITE(18,*) '#_Tpts indicates the number of points in the
* traced (_T) file'
WRITE(18,*) '#_Acpts indicates the number of points in the
* actual input file'
WRITE(18,*) 'Species is a code,
* #Ch is the # of chambers in
* last whorl, and coil is the
* coiling direction.'
WRITE(18,*) 'R (=X/Y) is the flatness
* of shell in keel view.'
WRITE(18,*) 'Rupkeel and Rlokeel are the absolute values of
* "upper and lower Kruemmungsradii" (in µm), respectively,
* of the ',TINTNORM
WRITE(18,*) ' '
WRITE(18,*) 'Sample, Age (Ma) #_Tpts #_Acpts
* Species #Ch Coil X (µm) Y (µm) R=X/Y Ar(mm2)
* Rupkeel Rlokeel'
C
END IF
C
C
C
C
WRITE(9,*) 'Write keel region coordinates to file ?'
WRITE(9,*) ' (yes = 1, no = 0)'
READ(9,*) ANSWER
C
C
C Recording errors during program run (Missing_files records the file
C names, that are listed in MEASUREMENTS, but that were not found during
C program run. The program continues then with the next file):
C
OPEN(19,FILE='Missing_files',STATUS='NEW')
C
C
C Input file handling:
C First get filenames of individual outlines from MEASUREMENTS, and the
C individual variables N through Area.
C MEASUREMENTS contains names of files with traced outlines
C (_T files, each name 17 chars long):
C
K=0
C
C*** (Main loop starts here):
C
C 8 READ(27,10,END=999) SPECIMEN,AGE,N,DX,DY,AREA ! 27 is file MEASUREMENTS
8 READ(27,*,END=999) SPECIMEN,AGE,N,DX,DY,AREA ! 27 is file MEASUREMENTS
C 10 FORMAT(A15,X,F5.2,X,I4,X,F6.1,X,F6.1,X,F7.4)
C
IF (TYPE.EQ.1) THEN
INPUT=SPECIMEN//'_T'
ELSE IF (TYPE.EQ.2) THEN
INPUT=SPECIMEN//'_INT'
ELSE IF (TYPE.EQ.3) THEN
INPUT=SPECIMEN//'_NORM'
END IF
C
WRITE(9,15) INPUT
C
15 FORMAT('. . .Sample ',A20,'. . .')
C
C
C Check, whether file with name INPUT exists or not. If it does not exist:
C print a message and continue reading next image.
C
INQUIRE(FILE=INPUT,EXIST=FILEX)
IF (.NOT.FILEX) THEN
WRITE(9,*) 'Error: file ',INPUT,' is missing'
WRITE(19,97) INPUT
97 FORMAT('Error: file ',A20,' is missing')
WRITE(9,*) '. . .Program continues. . .'
GOTO 8
END IF
C
C
C K counts the number of processed outlines (=control if all outlines
C are treated)
C
K=K+1
C
C
C Count the number of points in file INPUT (for output to file RESULT):
C
CALL COUNT(INPUT,Q)
C
C Next, find X,Y coordinates of upper and lower keel region
C by subroutines MAXIMA and MINIMA, respectively.
C Results are written to two external files filename_MaxK and filename_MinK:
C Top and base points of keel region are stored in X1,YMAX, and X2,YMIN, respectively.
C
CALL MAXIMA(E,INPUT,NMX,X1,YMAX)
CALL MINIMA(E,INPUT,NMI,X2,YMIN)
C
C
C X1 is the x-coordinate of the upper keel maximum point, where the x coord is the mean of all
C x values having the absolute Ymax.
C X2 is the x-coordinate of the lower keel minimum point, where the x coord is the mean of all
C x values having the absolute Ymin.
C
C
C Next, calculate third or fifth order polynom regression to approximate the
C keel region with subroutine CREG3 or CREG5
C The results for the upper keel region is stored in vector A,
C and the results for the lower keel region is stored in vector B.
C Vectors A and B have both 6 elements.
C
C INPUTMAX IS THE FILENAME HAVING THE upper keel region and INPUTMIN is the
C filename with the lower keel region. First calculate names of INPUTMAX and
C INPUTMIN from the character variable INPUT:
C
INPUTMAX=INPUT(1:15)//'_MaxK'
INPUTMIN=INPUT(1:15)//'_MinK'
C
C Now polynom regression for the upper and lower keel region;
C output of calculations depending on whether the output file is
C the simple one (RESULT) or the extended one (RES_...):
C
C***main IF statement:
C
IF (ANS.EQ.'1') THEN ! Simple output to file RESULT
C
IF (DEGREE.EQ.3) THEN
CALL CREG3(INPUTMAX,A,ERR) ! upper keel region
IF (ERR.EQ.'TFPIC') THEN
WRITE(18,78) INPUT,AGE,N,DX,DY,AREA
78 FORMAT(A20,', ',F5.2,', ',I4,', ',F6.1,', ',F6.1,', ',F7.4,
1' no Kruemmungskreis: too few points in keel region')
GOTO 95 ! advance to next sample.
END IF
CALL CREG3(INPUTMIN,B,ERR) ! lower keel region
ELSE IF (DEGREE.EQ.5) THEN
CALL CREG5(INPUTMAX,A,ERR) ! upper keel region
IF (ERR.EQ.'TFPIC') THEN ! TFPIC means Too Few Points In Curve
WRITE(18,78) INPUT,AGE,N,DX,DY,AREA
GOTO 95 ! advance to next sample.
END IF
CALL CREG5(INPUTMIN,B,ERR) ! lower keel region
END IF
C
C
C Interrupt this sample if ERR has the value 'ERROR' (this is the case if
C the determinant in SLE is zero. Write a message to RESULTS and continue with
C next sample:
C
IF (ERR.EQ.'ERROR') THEN
WRITE(18,81) INPUT,AGE,N,Q,DX,DY,AREA
81 FORMAT(A20,', ',F5.2,', ',I4,', ',I4,', ',F6.1,
1', ',F6.1,', ',F7.4,
2' no Kruemmungskreis: Determinant in SLE is zero')
C
GOTO 95 ! advance to next sample.
END IF
C
C
C Next, calculate bending at X1 and X2 with
C subroutine BEND3(V,X,RHO) or BEND5(V,X,RHO). BEND3 applies for a cubic spline regression
C and BEND5 applies to a fifth order polynom regression.
C V is the vector containing the coefficients of the third or fifth order polynomial
C regression line through the points of the keel region, X is the X coordinate where the
C keel is located, RHO is the Kruemmungsradius (here called bending) of the keel at X.
C The results are written to the external file RESULTS
C
IF (DEGREE.EQ.3) THEN
CALL BEND3(A,X1,RHOMAX) ! Upper keel region, cubic spline
CALL BEND3(B,X2,RHOMIN) ! Lower keel region, cubic spline
ELSE IF (DEGREE.EQ.5) THEN
CALL BEND5(A,X1,RHOMAX) ! Upper keel region
CALL BEND5(B,X2,RHOMIN) ! Lower keel region
END IF
C
C
C Output of all data from MEASUREMENTS and the results to file RESULT:
C RHOMAX,RHOMIN are the Kruemmungsradius in the upper and lower keel region, respectively.
C KMAX,KMIN are the Kruemmung (absolute value) in the upper and lower keel region, respectively
C
KMAX=ABS(1/RHOMAX)
KMIN=ABS(1/RHOMIN)
C
C Calculate y-coordinates of the centers of the upper (Y1) and lower (Y2)
C Kruemmungskreise
C
Y1=YMAX-ABS(RHOMAX)
Y2=YMIN+ABS(RHOMIN)
C
WRITE(18,20) INPUT,AGE,N,Q,DX,DY,AREA,ABS(RHOMAX),ABS(RHOMIN),
1KMAX,KMIN,X1,Y1,X2,Y2
20 FORMAT(A20,', ',F5.2,', ',I4,', ',I4,', ',
1F6.1,', ',F6.1,', ',F7.4,
2', ',F8.3,', ',F8.3,', ',F8.5,', ',F8.5,', ',F8.3,', ',F8.3,
3', ',F8.3,', ',F8.3)
C
C
C
C Keep or delete files INPUTMAX or INPUTMIN, that contain x,y coordinates
C of the keel region
C
95 IF (ANSWER.EQ.0) THEN
OPEN(25,FILE=INPUTMAX,STATUS='OLD')
CLOSE(25,STATUS='DELETE')
OPEN(26,FILE=INPUTMIN,STATUS='OLD')
CLOSE(26,STATUS='DELETE')
END IF
C
C
C
C Read next sample:
C
GOTO 8
C
C
C*** Main If statement:
C
ELSE IF (ANS.EQ.'2') THEN ! (extended output to file RES_...)
C
C First, calculate flatness R of shell: R=DX/DY
C
R=DX/DY
C
IF (DEGREE.EQ.3) THEN
CALL CREG3(INPUTMAX,A,ERR) ! upper keel region
IF (ERR.EQ.'TFPIC') THEN
WRITE(18,96) INPUT,AGE,N,SPECIES,NCHAMB,COIL,DX,DY,R,AREA
96 FORMAT(A20,', ',F5.2,', ',I4,', ',
1A7,', ',A1,', ',A1,', ',F6.1,', ',F6.1,', ',
2F6.4,', ',F7.4,
3' no Kruemmungskreis: too few points in keel region')
GOTO 99 ! advance to next sample.
END IF
CALL CREG3(INPUTMIN,B,ERR) ! lower keel region
ELSE IF (DEGREE.EQ.5) THEN
CALL CREG5(INPUTMAX,A,ERR) ! upper keel region
IF (ERR.EQ.'TFPIC') THEN ! TFPIC means Too Few Points In Curve
WRITE(18,96) INPUT,AGE,N,SPECIES,NCHAMB,COIL,DX,DY,R,AREA
GOTO 99 ! advance to next sample.
END IF
CALL CREG5(INPUTMIN,B,ERR) ! lower keel region
END IF
C
C
C Interrupt this sample if ERR has the value 'ERROR' (this is the case if
C the determinant in SLE is zero. Write a message to RESULTS and continue with
C next sample:
C
IF (ERR.EQ.'ERROR') THEN ! if the number of points (=Q) in INPUT is too low for polynom
C ! regression
WRITE(18,94) INPUT,AGE,N,Q,SPECIES,NCHAMB,COIL,DX,DY,R,AREA
94 FORMAT(A20,', ',F5.2,', ',I4,', ',I4,', ',
1A7,', ',A1,', ',A1,', ',F6.1,', ',F6.1,', ',
2F6.4,', ',F7.4,
3' no Kruemmungskreis: Determinant in SLE is zero')
C
GOTO 99 ! advance to next sample.
END IF
C
C
C
C Next, calculate bending at X1 and X2 with
C subroutine BEND3(V,X,RHO) or BEND5(V,X,RHO). BEND3 applies for a cubic spline regression
C and BEND5 applies to a fifth order polynom regression.
C V is the vector containing the coefficients of the third or fifth order polynomial
C regression line through the points of the keel region, X is the X coordinate where the
C keel is located, RHO is the Kruemmungsradius (here called bending) of the keel at X.
C The results are written to the external file RESULTS
C
IF (DEGREE.EQ.3) THEN
CALL BEND3(A,X1,RHOMAX) ! Upper keel region, cubic spline
CALL BEND3(B,X2,RHOMIN) ! Lower keel region, cubic spline
ELSE IF (DEGREE.EQ.5) THEN
CALL BEND5(A,X1,RHOMAX) ! Upper keel region
CALL BEND5(B,X2,RHOMIN) ! Lower keel region
END IF
C
C
C Output of all data from MEASUREMENTS and the results to file RESULT:
C RHOMAX,RHOMIN are the Kruemmungsradius in the upper and lower keel region, respectively.
C KMAX,KMIN are the Kruemmung (absolute value) in the upper and lower keel region, respectively
C
KMAX=ABS(1/RHOMAX)
KMIN=ABS(1/RHOMIN)
C
C Calculate y-coordinates of the centers of the upper (Y1) and lower (Y2)
C Kruemmungskreise
C
Y1=YMAX-ABS(RHOMAX)
Y2=YMIN+ABS(RHOMIN)
C
WRITE(18,30) INPUT,AGE,N,Q,SPECIES,NCHAMB,COIL,
1DX,DY,R,AREA,ABS(RHOMAX),ABS(RHOMIN)
30 FORMAT(A20,', ',F5.2,', ',I4,', ',I4,', ',
1A7,', ',A1,', ',A1,', ',
2F6.1,', ',F6.1,', ',F6.4,', ',F7.4,
3', ',F8.3,', ',F8.3)
C
C
C
C Keep or delete files INPUTMAX or INPUTMIN, that contain x,y coordinates
C of the keel region
C
99 IF (ANSWER.EQ.0) THEN
OPEN(25,FILE=INPUTMAX,STATUS='OLD')
CLOSE(25,STATUS='DELETE')
OPEN(26,FILE=INPUTMIN,STATUS='OLD')
CLOSE(26,STATUS='DELETE')
END IF
C
C
C
C Read next sample:
C
GOTO 8
C
END IF ! (End simple or extended output)
C
999 WRITE(9,*) K,' outlines processed'
PAUSE 1000
1000 CONTINUE
C
STOP
END

SUBROUTINE COUNT(INPUT,P)
C
C Determines the number of x,y coordinates (number of points) in file INPUT.
C P is returned to the calling unit.
C
INTEGER P
DOUBLE PRECISION X,Y
CHARACTER*20 INPUT
C
C Initialize P
C
P=0
C
OPEN(17,FILE=INPUT,STATUS='OLD')
5 READ(17,*,END=99) X,Y
P=P+1
GOTO 5
C
99 CLOSE(17)
C
RETURN
END

SUBROUTINE MAXIMA(EE,SAMPLE,NNMX,MXYMAX,MYMAX)
C
INTEGER I,EE,NNMX
DOUBLE PRECISION YMAX(50),XYMAX(50)
DOUBLE PRECISION X,Y,SUMX,MXYMAX,MYMAX
CHARACTER*20 SAMPLE
CHARACTER*20 OUTPUT
C
C Initialize XYMAX and YMAX:
C
DO 5, I=1,20
YMAX(I)=0.0
XYMAX(I)=0.0
5 CONTINUE
C
C Determine value of OUTPUT from variable SAMPLE:
C
OUTPUT=SAMPLE(1:15)//'_MaxK'
C
C Open Input and Output files:
C Input file is a file with name SAMPLE containing the cartesian X,Y coordinates of the curve
C Output file is MAXIMA containing the upper keel region X,Y coordinates.
C
OPEN(15,FILE=SAMPLE,STATUS='OLD')
OPEN(16,FILE=OUTPUT,STATUS='NEW')
C
C Calculating YMAX:
C (The X coordinates of the Y extrema points are stored in XYMAX).
C
C First find absolute maximum YMAX(1):
C
C WRITE(9,*) '. . . finding maxima . . .'
READ(15,*) X,YMAX(1)
10 READ(15,*,END=100) X,Y
IF (Y.GT.YMAX(1)) THEN
YMAX(1)=Y
XYMAX(1)=X
END IF
GOTO 10
C
100 REWIND(15)
C
C Check, whether there are more points with the same Y coordinate and
C write results to Output file MAXIMA:
C (The first value of YMAX(1), that was found from above is also printed here).
C
C Calculation of the X,Y coordinate of the location of the mean absolute
C maximum (mean of X[Y(1)],Y(1) ). MXYMAX denotes the mean X coordinate
C having the absolute Maximum.
C
C NNMX denotes the actual number of points, that belong to a particular
C maximum or submaximum.
C
SUMX=0.0
NNMX=0
12 READ(15,*,END=150) X,Y
IF (Y.EQ.YMAX(1)) THEN
NNMX=NNMX+1
SUMX=SUMX + X
WRITE(16,50) X,Y
50 FORMAT(F8.2,',',F8.2)
C
C Stop criterium: if there are more points exported than EE: Stop
C
IF (NNMX.EQ.EE) THEN
GOTO 300
END IF
END IF
GOTO 12
C
150 MXYMAX=SUMX/NNMX
MYMAX=YMAX(1)
C
REWIND(15)
C
C
C
C Now, finding subsequent submaxima:
C
DO 20, I=2,EE
READ(15,*,END=300) X,YMAX(I)
30 READ(15,*,END=200) X,Y
IF ((Y.GE.YMAX(I)).AND.(Y.LT.YMAX(I-1))) THEN
XYMAX(I)=X
YMAX(I)=Y
END IF
GOTO 30
C
200 REWIND(15)
C
C Again, check whether there are more points with the same Y coordinate and
C write results to output file MAXIMA:
C
14 READ(15,*,END=250) X,Y
IF (Y.EQ.YMAX(I)) THEN
NNMX=NNMX+1
WRITE(16,50) X,Y
C
C Stop criterium: if there are more points read than EE: Stop
C
IF (NNMX.EQ.EE) THEN
GOTO 300
END IF
END IF
GOTO 14
C
250 REWIND(15)
20 CONTINUE
C
300 CLOSE(15)
CLOSE(16)
C
RETURN
END

SUBROUTINE MINIMA(EE,SAMPLE,NNMI,MXYMIN,MYMIN)
C
INTEGER I,EE,NNMI
DOUBLE PRECISION YMIN(50),XYMIN(50)
DOUBLE PRECISION X,Y,SUMX,MXYMIN,MYMIN
CHARACTER*20 SAMPLE
CHARACTER*20 OUTPUT
C
C Initialize XYMIN and YMIN:
C
DO 5, I=1,20
YMIN(I)=0.0
XYMIN(I)=0.0
5 CONTINUE
C
C Determine value of OUTPUT from variable SAMPLE:
C
OUTPUT=SAMPLE(1:15)//'_MinK'
C
C Open Input and Output files:
C Input file is a file with name SAMPLE containing the cartesian X,Y coordinates of the curve
C Output file is MINIMA containing the lower keel region X,Y coordinates.
C
OPEN(15,FILE=SAMPLE,STATUS='OLD')
OPEN(17,FILE=OUTPUT,STATUS='NEW')
C
C YMIN:
C (The X coordinates of the Y extrema points are stored in XYMIN).
C
C First find absolute minimum YMIN(1):
C
C WRITE(9,*) '. . . finding minima . . .'
READ(15,*) X,YMIN(1)
310 READ(15,*,END=400) X,Y
IF (Y.LT.YMIN(1)) THEN
YMIN(1)=Y
XYMIN(1)=X
END IF
GOTO 310
C
400 REWIND(15)
C
C
C
C Check, whether there are more than one points with the same Y coordinate and
C write results to Output file MINIMA:
C (The first value of YMIN(1), that was found from above is also printed there).
C
C Calculation of the X,Y coordinate of the location of the mean absolute
C minimum (mean of X[Y(1)],Y(1) ). MXYMIN denotes the mean X coordinate
C having the absolute Minimum.
C
C NNMI denotes the actual number of points, that belong to a particular
C minimum or subminimum.
C
SUMX=0.0
NNMI=0
312 READ(15,*,END=450) X,Y
IF (Y.EQ.YMIN(1)) THEN
NNMI=NNMI+1
SUMX=SUMX + X
WRITE(17,50) X,Y
50 FORMAT(F8.2,',',F8.2)
C
C Stop criterium: if there are more points exported than EE: Stop
C
IF (NNMI.EQ.EE) THEN
GOTO 800
END IF
END IF
GOTO 312
C
450 MXYMIN=SUMX/NNMI
MYMIN=YMIN(1)
C
REWIND(15)
C
C
C
C
C Now, finding subsequent subminima:
C
DO 500, I=2,EE
READ(15,*,END=800) X,YMIN(I)
330 READ(15,*,END=600) X,Y
IF ((Y.LE.YMIN(I)).AND.(Y.GT.YMIN(I-1))) THEN
XYMIN(I)=X
YMIN(I)=Y
END IF
GOTO 330
C
600 REWIND(15)
C
C Again, check whether there are more points with the same Y coordinate and
C write results to output file MINIMA:
C
314 READ(15,*,END=550) X,Y
IF (Y.EQ.YMIN(I)) THEN
NNMI=NNMI+1
WRITE(17,50) X,Y
C
C Stop criterium: if there are more points read than EE: Stop
C
IF (NNMI.EQ.EE) THEN
GOTO 800
END IF
END IF
GOTO 314
C
550 REWIND(15)
500 CONTINUE
C
C
800 CLOSE(15)
CLOSE(17)
C
RETURN
END

SUBROUTINE CREG3(INPUT,V,ERR)
C
C Third order polynomial regression (cubic spline)
C
C Version 1.0, 1 January, 2001, by Michael Knappertsbusch
C Modified and adapted for implementation as a subroutine from program CURVREG5 version 1.0.
C
C INPUT:
C Given is a file with the name of character variable INPUT, which contains a suite of
C P points with cartesian coordinates Xi,Yi.
C The value of character variable INPUT enters the subroutine from the main unit.
C
C NOTE: For application of keel region approximation in globorotalid foraminifera:
C Input files should have a name of the form filename_MaxK or filename_MinK
C with 20 chars in order to follow the common file handling convention.
C
C The program determines a third order polynom approximation through these
C points by curvilinear regression.
C
C Condition: Number P of X,Y pairs > 4
C
C The regression line follows the equation
C
C y = aX**3 + bX**2 + cX + d
C
C The coefficients a through d can be determined by solving a system of
C 4 linear equations of the form
C
C U * X = V
C
C with X being the solution vector containing the coefficients a through d.
C Matrix U is defined as follows:
C
C Sum Xi**6 Sum Xi**5 Sum Xi**4 Sum Xi**3
C Sum Xi**5 Sum Xi**4 Sum Xi**3 Sum Xi**2
C Sum Xi**4 Sum Xi**3 Sum Xi**2 Sum Xi
C Sum Xi**3 Sum Xi**2 Sum Xi P
C
C
C Vector V is defined as:
C
C
C Sum Xi**3 Yi
C Sum Xi**2 Yi
C Sum Xi Yi
C Sum Yi
C
C With Sum meaning the summation of the summand from i=1 to P, and P being
C the number of X,Y pairs.
C
C OUTPUT:
C The vector V contains the coefficients a through d.
C Vector V is returned to the calling unit.
C
C
C File handling convention:
C All file names for morphometric analyses follow the same name convention.
C
C filename = Name of sample, 15 chars long.
C filename_T = File of raw curce with X,Y coordinates; 17 chars long.
C filename_POL = File with polar coordinates; 19 chars long.
C filename_INT = File with (cartesian) interpolated coords, 19 chars long.
C filename_MaxK = File with coords of upper keel region, 20 chars long.
C filename_MinK = File with coords of lower keel region, 20 chars long.
C filename_NORM = File with normalized curve, 20 chars long.
C
INTEGER I,J,M
DOUBLE PRECISION X,Y
DIMENSION U(4,4),V(4)
CHARACTER*5 ERR ! ERR indicates if the Determinant in SLE is zero.
CHARACTER*20 INPUT
PARAMETER(M=4) ! M indicates the number of equations to be solved.
C
C
C INITIALIZING U:
DO 5, I=1,4
DO 6, J=1,4
U(I,J)=0.0
6 CONTINUE
5 CONTINUE
C
C INITIALIZE V:
DO 7, I=1,4
V(I)=0.0
7 CONTINUE
C
C Open the INPUT file, read X,Y, and determine Elements of matrix U
C from sums of Xi coordinates.
C
C First open input file INPUT:
C
OPEN(20,FILE=INPUT,STATUS='OLD')
C
C Read X,Y coordinates, calculate sums, and determine number
C of points (P, is calculated in U(4,4):
C
10 READ(20,*,END=999) X,Y
C
C Components of matrix U: ! All other components of U remain zero.
C
U(1,1)=U(1,1) + X**6
U(1,2)=U(1,2) + X**5
U(1,3)=U(1,3) + X**4
U(1,4)=U(1,4) + X**3
U(2,1)=U(2,1) + X**5
U(2,2)=U(2,2) + X**4
U(2,3)=U(2,3) + X**3
U(2,4)=U(2,4) + X**2
U(3,1)=U(3,1) + X**4
U(3,2)=U(3,2) + X**3
U(3,3)=U(3,3) + X**2
U(3,4)=U(3,4) + X
U(4,1)=U(4,1) + X**3
U(4,2)=U(4,2) + X**2
U(4,3)=U(4,3) + X
U(4,4)=U(4,4) + 1.0

C
C Components of vector V: ! All other components of V remain zero.
C
V(1)=V(1) + Y * X**3
V(2)=V(2) + Y * X**2
V(3)=V(3) + Y * X
V(4)=V(4) + Y
C
GOTO 10
C
999 CONTINUE
C
C
C Stop criterion if the condition (number of points > 4) is not fulfilled:
C
IF (U(4,4).LT.4.0) THEN
WRITE(9,*) 'Too few points in curve.
* Program advances to next sample'
ERR='TFPIC' ! TFPIC means: Too Few Points In Curve
RETURN
END IF
C
CLOSE(20)
C
C
C Now solve system with 4 equations. Use of program code from SEQ.f
C and subroutine SLE, which were adapted to N=4.
C The results are returned in vecor V:
C
CALL SLE(U,V,1.0E-07,M,ERR)
C
RETURN
END

SUBROUTINE CREG5(INPUT,V,ERR)
C
C Fifth order polynomial regression
C
C Version 2.0, 21 August, 2000, by Michael Knappertsbusch
C Modified and adapted for implementation as a subroutine from program CURVREG5 version 1.0.
C
C INPUT:
C Given is a file with the name of character variable INPUT, which contains a suite of
C P points with cartesian coordinates Xi,Yi.
C The value of character variable INPUT enters the subroutine from the main unit.
C
C NOTE: For application of keel region approximation in globorotalid foraminifera:
C Input files should have a name of the form filename_MaxK or filename_MinK
C with 20 chars in order to follow the common file handling convention.
C
C The program determines a fifth order polynom approximation through these
C points by curvilinear regression.
C
C Condition: Number P of X,Y pairs > 6
C
C The regression line follows the equation
C
C y = aX**5 + bX**4 + cX**3 + dX**2 + eX + f
C
C The coefficients athrough f can be determined by solving a system of
C 6 linear equations of the form
C
C U * X = V
C
C with X being the solution vector containing the coefficients a through f.
C Matrix U is defined as follows:
C
C Sum Xi**10 Sum Xi**9 Sum Xi**8 Sum Xi**7 Sum Xi**6 Sum Xi**5
C Sum Xi**9 Sum Xi**8 Sum Xi**7 Sum Xi**6 Sum Xi**5 Sum Xi**4
C Sum Xi**8 Sum Xi**7 Sum Xi**6 Sum Xi**5 Sum Xi**4 Sum Xi**3
C Sum Xi**7 Sum Xi**6 Sum Xi**5 Sum Xi**4 Sum Xi**3 Sum Xi**2
C Sum Xi**6 Sum Xi**5 Sum Xi**4 Sum Xi**3 Sum Xi**2 Sum Xi
C Sum Xi**5 Sum Xi**4 Sum Xi**3 Sum Xi**2 Sum Xi P
C
C
C Vector V is defined as:
C
C Sum Xi**5 Yi
C Sum Xi**4 Yi
C Sum Xi**3 Yi
C Sum Xi**2 Yi
C Sum Xi Yi
C Sum Yi
C
C With Sum meaning the summation of the summand from i=1 to P, and P being
C the number of X,Y pairs.
C
C OUTPUT:
C The vector V contains the coefficients a through f.
C Vector V is returned to the calling unit.
C
C
C File handling convention:
C All file names for morphometric analyses follow the same name convention.
C
C filename = Name of sample, 15 chars long.
C filename_T = File of raw curce with X,Y coordinates; 17 chars long.
C filename_POL = File with polar coordinates; 19 chars long.
C filename_INT = File with (cartesian) interpolated coords, 19 chars long.
C filename_MaxK = File with coords of upper keel region, 20 chars long.
C filename_MinK = File with coords of lower keel region, 20 chars long.
C filename_NORM = File with normalized curve, 20 chars long.
C
INTEGER I,J,M
DOUBLE PRECISION X,Y
DIMENSION U(6,6),V(6)
CHARACTER*5 ERR ! ERR indicates if the Determinant in SLE is zero or if there are
C ! too many points.
CHARACTER*20 INPUT
PARAMETER(M=6) ! M indicates the number of equations to be solved.
C
C
C!!! WRITE(9,*) '. . .Program determines coefficients for. . .'
C!!! WRITE(9,*) ' a fifth order polynomial regression'
C!!! WRITE(9,*) ' through a set of points'
C!!! WRITE(9,*) ' '
C!!! WRITE(9,*) ' (single user mode)'
C!!! WRITE(9,*) 'The number of points available must be > 6'
C
C
C INITIALIZING U:
DO 5, I=1,6
DO 6, J=1,6
U(I,J)=0.0
6 CONTINUE
5 CONTINUE
C
C INITIALIZE V:
DO 7, I=1,6
V(I)=0.0
7 CONTINUE
C
C Open the INPUT file, read X,Y, and determine Elements of matrix U
C from sums of Xi coordinates.
C
C First open input file INPUT:
C
OPEN(20,FILE=INPUT,STATUS='OLD')
C
C Read X,Y coordinates, calculate sums, and determine number
C of points (P, is calculated in U(6,6):
C
10 READ(20,*,END=999) X,Y
C
C Components of matrix U:
C
U(1,1)=U(1,1) + X**10
U(1,2)=U(1,2) + X**9
U(1,3)=U(1,3) + X**8
U(1,4)=U(1,4) + X**7
U(1,5)=U(1,5) + X**6
U(1,6)=U(1,6) + X**5
U(2,1)=U(2,1) + X**9
U(2,2)=U(2,2) + X**8
U(2,3)=U(2,3) + X**7
U(2,4)=U(2,4) + X**6
U(2,5)=U(2,5) + X**5
U(2,6)=U(2,6) + X**4
U(3,1)=U(3,1) + X**8
U(3,2)=U(3,2) + X**7
U(3,3)=U(3,3) + X**6
U(3,4)=U(3,4) + X**5
U(3,5)=U(3,5) + X**4
U(3,6)=U(3,6) + X**3
U(4,1)=U(4,1) + X**7
U(4,2)=U(4,2) + X**6
U(4,3)=U(4,3) + X**5
U(4,4)=U(4,4) + X**4
U(4,5)=U(4,5) + X**3
U(4,6)=U(4,6) + X**2
U(5,1)=U(5,1) + X**6
U(5,2)=U(5,2) + X**5
U(5,3)=U(5,3) + X**4
U(5,4)=U(5,4) + X**3
U(5,5)=U(5,5) + X**2
U(5,6)=U(5,6) + X
U(6,1)=U(6,1) + X**5
U(6,2)=U(6,2) + X**4
U(6,3)=U(6,3) + X**3
U(6,4)=U(6,4) + X**2
U(6,5)=U(6,5) + X
U(6,6)=U(6,6) + 1.0
C
C Components of vector V:
C
V(1)=V(1) + Y * X**5
V(2)=V(2) + Y * X**4
V(3)=V(3) + Y * X**3
V(4)=V(4) + Y * X**2
V(5)=V(5) + Y * X
V(6)=V(6) + Y
C
GOTO 10
C
999 CONTINUE
C
C Stop criterion if the condition (number of points > 6) is not fulfilled:
C
IF (U(6,6).LT.6.0) THEN
WRITE(9,*) 'Too few points in curve.
* Program advances to next sample'
ERR='TFPIC' !TFPIC means: Too Few Points In Curve
RETURN
END IF
C
CLOSE(20)
C
C
C Now solve system with 6 equations. Use of program code from SEQ.f
C and subroutine SLE.
C The results are returned in vecor V:
C
CALL SLE(U,V,1.0E-07,M,ERR)
C
C
RETURN
END

SUBROUTINE SLE(A,B,ZERO,N,ERR)
C
C SUBROUTINE FOR SOLUTION OF 6 SIMULTANEOUS EQUATIONS.
C MATRIX A IS N X N AND B IS A COLUMN VECTOR OF N ELEMENTS.
C Adapted from J.C. Davis (1973) Statistics and Data Analysis
C in Geology.
C
C EQUATION OF THE FORM A * X = B
C
C A CONTAINS THE COEFFICIENTS, B CONTAINS THE CONSTANTS,
C X CONTAINS UNKNOWN VARIABLE X.
C
C A IS CONVERTED TO THE IDENTITY MATRIX.
C VECTOR B ENTERS THE SUBROUTINE WITH THE CONSTANTS, AND
C VECTOR B IS RETURNED TO THE CALLING UNIT CONTAINING THE
C VARIABLE X (THIS MEANS THAT B RETURNS THE SOLUTION).
C
INTEGER N ! N is the number of equations to be solved (variable dimensioning).
DOUBLE PRECISION ZERO,DIV,RATIO
CHARACTER*5 ERR ! ERR indicates in the Determinant is zero.
DIMENSION A(N,N),B(N)
C
C
DO 100 I=1,N
DIV=A(I,I)
IF (ABS(DIV)-ZERO) 99,99,1
1 DO 101 J=1,N
A(I,J)=A(I,J)/DIV
101 CONTINUE
C
B(I)=B(I)/DIV
DO 102 J=1,N
IF (I-J) 2,102,2
2 RATIO=A(J,I)
DO 103 K=1,N
A(J,K)=A(J,K)-RATIO*A(I,K)
103 CONTINUE
C
B(J)=B(J)-RATIO*B(I)
102 CONTINUE
100 CONTINUE
C
RETURN
C
99 WRITE(9,*) 'Determinant of matrix is zero: no single solution.'
WRITE(9,*) '(-->Either no solution or indefinitely many
* solutions. Program reads next sample.)'
ERR='ERROR'
C
RETURN
END

SUBROUTINE BEND3(V,X,BEND)
C
C Subroutine to calculate the bending of a 3rd order polynomial at location X.
C This subroutine is applied to calculate the curvature of the outline of a
C globorotalid foraminifera in the upper and lower keel region.
C
C By Michael Knappertsbusch, 4 January 2001
C Version 1.0
C
C
C Input is a vector with 6 components containing the coefficients of the 3rd order
C polynomial regression curve of the form
C
C y = aX**3 + bX**2 + cX + d
C
C and the variable X, which is located in the middle of the upper or
C lower keel region. The coefficients a to d are stored in vector V as
C elements V(1) to V(4), respectively (elements V(5) and V(6) are zero).
C
C Output is the double precision variable BEND (=bending, or curvature at X), which is
C returned to the calling unit.
C
C
C The bending is calculated by the formula:
C
C [1 + (f'(X))**2]**3/2
C BEND(x) = ----------------------
C f''(X)
C
C f'(x) = 3*a*X**2 + 2*b*X + c
C f''(x) = 6*a*x + 2*b
C
DOUBLE PRECISION X,BEND
DOUBLE PRECISION SUM1,SUM2,SUM3
DIMENSION V(6)
C
C
C Calculation of partial sums to calculate the bending:
C
C
SUM1=3.0*V(1)*X**2 + 2.0*V(2)*X + V(3)
SUM2=6.0*V(1)*X + 2.0*V(2)
SUM3=1.0 + SUM1**2
BEND=(SUM3**3/2)/SUM2
C
RETURN
END

SUBROUTINE BEND5(V,X,BEND)
C
C Subroutine to calculate the bending of a 5th order polynomial at location X.
C This subroutine is applied to calculate the curvature of the outline of a
C globorotalid foraminifera in the upper and lower keel region.
C
C By Michael Knappertsbusch, 21. August 2000
C Version 1.0
C
C
C Input is a vector with 6 components containing the coefficients of the 5th order
C polynomial regression curve of the form
C
C y = aX**5 + bX**4 + cX**3 + dX**2 + eX + f
C
C and the variable X, which is located in the middle of the upper or
C lower keel region. The coefficients a to f are stored in vector V as
C elements V(1) to V(6), respectively.
C
C Output is the double precision variable BEND (=bending, or curvature at X), which is
C returned to the calling unit.
C
C

C The bending is calculated by the formula:
C
C [1 + (f'(X))**2]**3/2
C BEND(x) = ----------------------
C f''(X)
C
C
DOUBLE PRECISION X,BEND
DOUBLE PRECISION SUM1,SUM2,SUM3,SUM4
DOUBLE PRECISION SUM5,SUM6,SUM7,SUM8
DOUBLE PRECISION SUM9,SUM10,SUM11
DIMENSION V(6)
C
C
C Calculation of partial sums to calculate the bending:
C
SUM1=5*V(1)*X**4
SUM2=4*V(2)*X**3
SUM3=3*V(3)*X**2
SUM4=2*V(4)*X
SUM5=SUM1+SUM2+SUM3+SUM4+V(5)
SUM6=1+SUM5**2
C
SUM7=20*V(1)*X**3
SUM8=12*V(2)*X**2
SUM9=6*V(3)*X
SUM10=2*V(4)
SUM11=SUM7+SUM8+SUM9+SUM10
C
BEND=(SUM6**(3/2))/SUM11
C
RETURN
END

6.6. Listing of Program Norm52.f

PROGRAM NORM
C
C Version 5.3,
C Modified from version 5 from November 2, 1998, by Michael Knappertsbusch
C Modified from version 5.2
C
C Modifications:
C 19.12.2000: Addition of error reports for the cases if files are missing.
C 25.1.2004: Addition of Maximum Diameter Normalization.
C 10.8.2004: Some comments added.
C 13.9.2005: Option II in version 5.2 did a wrong normalization. This is now corrected for and
C some comments were added.
C
C Normalization of closed plane curves:
C Given are the cartesian coordinates of a closed plane curve of a particle silhouette,
C with the origin of the coordinate system being at the centroid of the curve (e.g. the
C output (_INT files) of program Sprep).
C
C In the program it is assumed, that the particles have all been previously oriented according to
C a common imaging protocol, so that the longest axis is vertical.
C
C Program Norm normalizes this curve using two methods, which can be
C selected by the user:
C
C Method I (=Mean Radius Normalization):
C The program converts the set of X and Y coordinates of the
C curve into polar coordinates R and Theta, calculates
C the AVERAGE radius RMEAN, and normalizes the individual values by division of R by RMEAN.
C Thereafter it reconverts them back into cartesian coordinates.
C !!!!! The output files with the mean radius normalized coordinates have the extension _NORM.
C The Mean Radius Normalization is used for Fourier analysis (where, for example, the zero'th
C harmonic function is a circle with radius 1).
C
C Method II (=Maximum Diameter Normalization):
C In Option II the program normalizes the x,y coordinates of the curve by the MAXIMUM diameter D
C of the curve. Because the specimens are all oriented according to the above give protocol, the
C longest axis of the particle is vertical. After normalization, this longest axis becomes equal
C to 1.
C !!!!! The output files with the normalized coordinates from option II have the extension _N in
C order to distinguish them from the output files from Method I.
C The Maximum Diameter Normalization is used for the geometric comparison of outline curves.
C
C
C !!! Processing in batch mode. !!!
C
C Input files are:
C - File with name FILE_LIST containing a list of all
C names of files (filename), that contain the X,Y coordinates
C to be normalized. FILE_LIST is a character variable
C and can bear any name of 20 characters length.
C
C - All files, that contain the X,Y coordinates of the
C curve to be normalized. The names of these files are
C interpreted through the character variable INPUT. The
C expression of INPUT is used to compose the name of the
C corresponding output file (suffix _NORM for option I and _N for option II).
C
C Output files are:
C - One single file for every input curve with the same name as the
C the input curve file (suffix _T), and carrying the suffix _NORM or _N.
C
C Note: All file handling for morphometric analyses follows a name convention:
C
C filename = Name of sample, 15 characters long.
C filename_T = Input file for raw curves with cartesian coordinates, 17 characters long.
C filename_POL = Polar coordinates of filename, 19 characters long.
C filename_INT = Interpolated (cartesian) coordinates of filename, 19 characters long.
C filename_NORM = Normalized curve, 20 characters long.
C
INTEGER N,NMAX
INTEGER ANSWER
LOGICAL FILEX
DOUBLE PRECISION X,Y,R,THETA,U,V
DOUBLE PRECISION SUM,RMEAN,NORMR
DOUBLE PRECISION YMAX,YMIN,D,XNORD,YNORD
CHARACTER*11 HEADER
CHARACTER*20 FILE_LIST,OUTPUT
CHARACTER*19 INPUT
C
WRITE(9,*) '************************************'
WRITE(9,*) '* *'
WRITE(9,*) '* Program Norm Version 5.3 *'
WRITE(9,*) '* *'
WRITE(9,*) '* by Michael Knappertsbusch *'
WRITE(9,*) '* *'
WRITE(9,*) '************************************'
WRITE(9,*) ' '
C WRITE(9,*) 'Preparation for outline analysis: Normalization:'
WRITE(9,*) ' '
WRITE(9,*) 'Note, that specimens must be vertically oriented'
WRITE(9,*) 'and that input files must have the centroid at'
WRITE(9,*) 'the origin of the coordinate system.'
WRITE(9,*) ' '
WRITE(9,*) 'Enter the method of normalization:'
WRITE(9,*) ' Mean Radius Normalization (1), or'
WRITE(9,*) ' Maximum Diameter Normalization (2)'
READ(9,*) ANSWER ! Selects for option I (=1) or II (=2)
WRITE(9,*) ' '
IF (ANSWER.EQ.1) THEN
WRITE(9,*) 'Mean Radius Normalization, output=_NORM'
ELSE IF (ANSWER.EQ.2) THEN
WRITE(9,*) 'Maximum Radius Normalization, output=_N'
END IF
C

C
C*** Record errors during program run:
C
OPEN(25,FILE='Errors_during_Norm',STATUS='NEW')
C
C
C*** Get filenames from file list FILE_LIST, open individual files
C and read X,Y data.
C
C WRITE(9,*) '. . .Names of interpolated outline data needed:. . .'
WRITE(9,*) 'enter list containing names of INT files (20 chars)'
READ(9,2) FILE_LIST
2 FORMAT(A20)
C
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_Norm":
C
READ(20,3) INPUT
HEADER=INPUT(1:11)
WRITE(25,7) HEADER
7 FORMAT('Sample name: ',A11)
REWIND(20)
C
C
C*** Main loop starts at statement #1.
C
C
1 READ(20,3,END=999) INPUT ! New curve is being read
WRITE(9,4) '. . .Reading new file ',INPUT,'. . .'
3 FORMAT(A19)
4 FORMAT(A22,A19,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
8 FORMAT('File ',A19,' is missing')
WRITE(9,*) '. . .Program continues. . .'
GOTO 1
END IF
C
C
C Determination of names for output files:
C
IF (ANSWER.EQ.1) THEN
OUTPUT=INPUT(1:15)//'_NORM' ! Output for Mean radius normalization
ELSE IF (ANSWER.EQ.2) THEN
OUTPUT=INPUT(1:15)//'_N' ! Output for maximum Diameter normalization
END IF
C
OPEN(15,FILE=INPUT,STATUS='OLD')
OPEN(16,FILE=OUTPUT,STATUS='NEW')
C
C
C
C Calculations:
C Mean Radius Normalization (Method I):
C
IF (ANSWER.EQ.1) THEN
C First calculating the mean radius RMEAN:
C
N=0
SUM=0.0
5 READ(15,*,END=10) X,Y
N=N+1
CALL POLAR(X,Y,R,THETA)
SUM=SUM+R
GOTO 5
C
10 NMAX=N
RMEAN=SUM/NMAX
C
WRITE(9,*) 'MEAN RADIUS IS: ',RMEAN
C
REWIND(15)
C
C
C Now calculating the normalized values of X and Y
C according to the selected method for standardization:
C
15 READ(15,*,END=20) X,Y
CALL POLAR(X,Y,R,THETA)
NORMR=R/RMEAN
C
C Reconvert into cartesian coord-system
C
U=NORMR*COS(THETA)
V=NORMR*SIN(THETA)
WRITE(16,*) U,',',V
GOTO 15
C
20 WRITE(9,*) INPUT,' normalized'
WRITE(9,*) 'File written with ',NMAX,' values'
CLOSE(15)
CLOSE(16)
END IF
C
C
C
C
C
C
C Maximum Radius normalization (Method II):
C
IF (ANSWER.EQ.2) THEN
N=0
C Determination of Ymax, Ymin, and then the vertical distance D=Ymax+Abs(Ymin):
C
READ(15,*) X,Y
N=N+1
YMAX=Y
30 READ(15,*,END=40) X,Y
N=N+1
IF (Y.GT.YMAX) THEN
YMAX=Y
END IF
GOTO 30
C
40 NMAX=N
REWIND(15)
C
C
READ(15,*) X,Y
YMIN=Y
45 READ(15,*,END=50) X,Y
IF (Y.LT.YMIN) THEN
YMIN=Y
END IF
GOTO 45
C
50 REWIND(15)
C
D=YMAX+ABS(YMIN)
C
C
C
C Now normalize each x,y pair with respect to D:
C
70 READ(15,*,END=80) X,Y
XNORD=X/D
YNORD=Y/D
WRITE(16,*) XNORD,',',YNORD
GOTO 70
C
80 WRITE(9,*) INPUT,' normalized'
WRITE(9,*) 'File written with ',NMAX,' values'
CLOSE(15)
CLOSE(16)
END IF
C
C
C
GOTO 1
C
999 WRITE(9,*) '. . .All files treated. . .'
CLOSE(25)
PAUSE 1000
1000 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=2*PI/3
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

6.7. Listing for program Compose31.f

PROGRAM COMPOSE
C
C Version 3.1 by Michael Knappertsbusch, 19 December 2003
C
C !!! batch processing mode !!!
C
C This program generates a data file with morphometric measurements
C from single specimens, that were created with program KeelWidth93.
C
C Input files:
C 1.) FILE_LIST, contains the names of all individual files to be treated.
C 2.) The individual files with the measurements of each specimen. The files,
C that contain the measurements have the general name
C "RESKW_sample_species_nChamb_COIL". They were generated with program
C KeelWidt93.out, and contain the following variables in this order:
C
C SPECIMEN,AGE,N,SPECIES,NCHAMB,COIL,DX,DY,R,AREA,D10,D90,D9010,DMIN,DMAX,PHI1,
C PHI2,ARESLO,ARESUP,DEV
C
C Output files:
C One file with the name "Composed_files" is generated.
C In each file "RESKW_sample_species_nChamb_COIL" the program strips off
C the header information and writes the morphometric
C data into "Composed_files". This file can be read with a commercial data
C plotting program (Cricket Graph, Excel) to generate scatter plots.
C
C SPECIMEN: = Name of the specimen in a sample
C AGE: = Age, in million years
C N: = Number of points in original outline
C SPECIES: = Name of the species (given as a 7 chars long code)
C NCHAMB: = Number of chambers in last whorl
C COIL: = Coiling direction
C DX,DY: = Width (spiral height) and length in keel view (in µm)
C R: = Ratio of DX/DY
C AREA: = Area in keel view, in mm2
C D10: = Lower relative keelwidth (at 10% axial extension), expressed as fraction of the width of
C the shell (D10/DX).
C D90 = Upper relative keelwidth (at 90% axial extension), expressed as fraction of the width of
C the shell (D90/DX).
C D9010 = Asolute difference of (D90-D10), expressed as fraction of the width of the shell
C [ABS(D90-D10)/DX].
C DMIN: = DMIN=MIN(D10,D90)/DX
C DMAX: = DMAX=MAX(D10,D90)/DX
C PHI1: = Upper keel angle, in degrees
C PHI2: = Lower keel angle, in degrees
C ARESLO: = Lower keel residual area, in mm2
C ARESUP: = Upper keel residual area, in mm2
C DEV: = Deviation of total keel residual area from the value of AREA, in per cent.
C
C
INTEGER N
REAL AGE,DX,DY,R,AREA,D10,D90
REAL D9010,DMIN,DMAX,PHI1,PHI2
REAL ARESLO,ARESUP,DEV
CHARACTER*1 COIL
CHARACTER*2 NCHAMB
CHARACTER*7 SPECIES
CHARACTER*17 SPECIMEN
CHARACTER*20 FILE_LIST ! FILE_LIST reads the name of the file containing
C the INPUT files.
CHARACTER*1 HEADER1
CHARACTER*31 INPUT ! Reads the filenames of "RESKW_sample_species_nChamb_COIL"
C
C
C
WRITE(9,*) '*****************************************'
WRITE(9,*) '* *'
WRITE(9,*) '* Program Compose (preparation for plotting data) *'
WRITE(9,*) '* . . .batch processing mode. . . *'
WRITE(9,*) '* *'
WRITE(9,*) '*Version 3.1 by Michael Knappertsbusch,16.12.2003 *'
WRITE(9,*) '* Only use in combination with KeelWidth93.out *'
WRITE(9,*) '* *'
WRITE(9,*) '******************************************'
C
C First, open output file, and write a header line in it so, that
C it can be read by Cricket Graph:
C
OPEN(14,FILE='Composed_files',STATUS='NEW')
WRITE(14,*) '*'
WRITE(14,2)
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 %)')
C
C
C Now, get name of the file list containing the input files,
C and start reading the INPUT file names.
C
WRITE(9,*) 'Enter the name of the list with
1 all samples (max 20 chars long)'
READ(9,'(A20)') FILE_LIST
C
C
C Open FILE_LIST, read an INPUT file name from it:
OPEN(13,FILE=FILE_LIST,STATUS='OLD')
C
C*** Main sample loop starts here:
C
5 READ(13,'(A31)',END=999) INPUT
WRITE(9,6) INPUT
6 FORMAT('. . .Sample ',A31,'. . .')
C
C Open the INPUT file, and read the header block (first 6 lines) in it:
C
OPEN(15,FILE=INPUT,STATUS='OLD')
C
DO 10 I=1,6
READ(15,'(A1)') HEADER1
10 CONTINUE
C
C Read subsequent data lines in that file:
C
20 READ(15,30,END=99) SPECIMEN,AGE,N,SPECIES,
1 NCHAMB,COIL,DX,DY,R,AREA,D10,D90,
2 D9010,DMIN,DMAX,PHI1,PHI2,
3 ARESLO,ARESUP,DEV
C
30 FORMAT(A17,2X,F5.2,2X,I4,2X,A7,
1 1X,A2,2X,A1,2X,F6.1,2X,F6.1,
2 2X,F6.4,2X,F7.4,2X,F7.2,2X,F7.2,
3 2X,F7.2,2X,F7.2,2X,F7.2,2X,F7.3,
4 2X,F7.3,2X,F7.4,7X,F7.4,
5 6X,F7.3)
C
C
C Write data to output file:
C
WRITE(14,40) SPECIMEN,AGE,N,SPECIES,
1 NCHAMB,COIL,DX,DY,R,AREA,D10,D90,
2 D9010,DMIN,DMAX,PHI1,PHI2,
3 ARESLO,ARESUP,DEV
40 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)
C
GOTO 20 ! Read next data line
C
99 CLOSE(15) ! Close actual sample file
GOTO 5 ! Read next sample in FILE_LIST
C
999 WRITE(9,*) 'One file "Composed_files" created'
PAUSE 1000
1000 CONTINUE
STOP
END

6.8. Listing for program Stat_Prep21.f

program Stat_Prep
C
C Version 2.1, 19.12.2003, by Michael Knappertsbusch
C modified from Version 1.0
C
C !!! Note: single processing mode !!!
C And: Programmed for the extended output option from KeelWidth93.out
C
C This program prepares the input file for program Stats61.out or Stats62.out if
C split-weighted averages need to be calculated with stats61.out (batch) during
C cycle II of the morphometric analysis.
C
C Split-weighted averaging is necessary if in a single sample foraminiferal specimens
C were unevenly distributed in the different size classes.
C
C Example: Only 2 specimens of G. menardii "A" occur in entire split of size class 500-1000µm,
C but 50 specimens of the same species occur in 1/32 split of size class 100-500µm.
C
C In this example, the statistical weight of the larger size fraction is very small, and
C average size parameters are dominated by the smaller specimens. Here, program Stat_prep
C rearranges the input data so that in the above example each of the two specimens from the
C size fraction 500-1000µm appear only as a single record (1/1 split), but each specimen of
C the 100-500µm size fraction appears 32 times in the output file. The output file is then
C directly fed to program Stats61.out (batch) or Stats62.out (batch).
C
C Input file:
C There is one single input file, usually file RESKW_SAMPLE_SPECIES_NCHAMB_COIL, which
C was generated with program KeelWidth93.out, and edited to add the split number before it
C is entered to program Stat_Prep21.out.
C
C Editing of RESKW_SAMPLE_SPECIES_NCHAMB_COIL can be done by hand in the following way:
C In the first column of the data body the split number must be inserted (=SPLIT), which
C is an integer indicating the factor by which each record must be pasted into the file
C (in the above example SPLIT is 1 for 1/1 splits, and 32 for 1/32 split).
C Next, a character variable SFRACT must be inserted (10 characters long) indicating from
C which size fraction the specimen comes from. Then, the databody follows, which is normally
C fed to Stats6.out (batch) if no split-weighted averaging would be necessary (e.g. specimen
C code, age, original number of points in outline, #Ch, Coil, etc.
C
C Example:
C FORMAT(I2,1X,A10,2X,A17,2X,F5.2,2X,I4,2X,A7,....)
C 32, 100-500µm, 502A211057K0602_NORM, 3.39, 585, menardB, 5, d, 268.8,......
C
C Output file:
C One output file is generated, that has precisely the same format as the original (unedited)
C version of file RESKW_SAMPLE_SPECIES_NCHAMB_COIL, but with multiple records of the samples according
C to the split number. The output file can directly be fed to program Stats6.out for further
C processing.
C
C
INTEGER SPLIT,TPTS ! SPLIT=Split factor, TPTS=# of points in original outline
REAL AGE,X(1:14) ! X(1) to X(14) are the morphometric parameters (X, Y, R, Ar,
C D10,D90,D9010,DMIN,DMAX,PHI1,PHI2,ARESLO,ARESUP, and DEV)
CHARACTER*10 SFRACT ! Describes the size fraction, from which a specimen comes from
CHARACTER*17 SPECIMEN ! Code for the specimens
CHARACTER*7 SPECIES ! Species names
CHARACTER*1 COIL ! COIL=coiling direction.
CHARACTER*2 NCHAMB ! NCHAMB=# Chambers
CHARACTER*31 INPUT,OUTPUT ! Variables for input and output files
CHARACTER*183 HEADER ! Reads the header info in the input file
C
C
C
WRITE(9,*) '**************************'
WRITE(9,*) '* *'
WRITE(9,*) '* Program Stat_Prep *'
WRITE(9,*) '* *'
WRITE(9,*) '* Version 2.1, 19.12.2003 *'
WRITE(9,*) '* By Michael Knappertsbusch *'
WRITE(9,*) '* *'
WRITE(9,*) '* !!! Single processing mode !!! *'
WRITE(9,*) '* *'
WRITE(9,*) '* Prepares the format of a file *'
WRITE(9,*) '* with split-weighted samples *'
WRITE(9,*) '* for input into program Stats *'
WRITE(9,*) '* *'
WRITE(9,*) '* Use only in combination with*'
WRITE(9,*) '* KeelWidth93.out, Stats61.out, *'
WRITE(9,*) '* or Stats62.out (=with Tabs) *'
WRITE(9,*) '* *'
WRITE(9,*) '**************************'
C
C
C
WRITE(9,*) ' '
WRITE(9,*) '. . .Example for data format. . .'
WRITE(9,*) '(I2,1X,A10,2X,A17,2X,F5.2,2X,I4,2X,A7,...)'
WRITE(9,*) '32, 100-500µm, 502A211057K0602_NORM, 3.39
*, 585, menardB, 5, d, 268.8'
WRITE(9,*) ' '
WRITE(9,*) 'Enter the name of the input file
1 (max 31 chars)'
READ(9,5) INPUT
5 FORMAT(A31)
C
C Determine name of output file:
C
WRITE(9,*) 'Enter the name of the output file
1 (max 31 chars)'
READ(9,5) OUTPUT
C
OPEN(15,FILE=INPUT,STATUS='OLD')
OPEN(16,FILE=OUTPUT,STATUS='NEW')
C
C Write same header information to output file as is given in the input file
C (there are 6 header lines that remain the same in the input and output files):
C
DO 10, I=1,6
READ(15,15) HEADER
WRITE(16,15) HEADER
10 CONTINUE
15 FORMAT(A183)
C
C
C
C Now, reading data lines from the INPUT file:
C
17 READ(15,20,ERR=998,END=999) SPLIT,SFRACT,SPECIMEN,AGE,
* TPTS,SPECIES,NCHAMB,COIL,X(1),X(2),
* X(3),X(4),X(5),X(6),X(7),X(8),X(9),X(10),
* X(11),X(12),X(13),X(14)
C
20 FORMAT(I2,1X,A10,2X,
* A17,2X,F5.2,2X,I4,2X,A7,1X,A2,2X,A1,2X,
* F6.1,2X,F6.1,2X,F6.4,2X,F7.4,2X,F7.2,2X,
* F7.2,2X,F7.2,2X,F7.2,2X,F7.2,2X,F7.3,2X,
* F7.3,2X,F7.4,7X,F7.4,6X,F7.3)
C
C ...and write to the OUTPUT file in the necessary format:
C
DO 30,I=1,SPLIT
WRITE(16,35) SPECIMEN,AGE,TPTS,SPECIES,
1NCHAMB,COIL,(X(I),I=1,14)
30 CONTINUE
C
35 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)
C
GOTO 17
C
998 WRITE(9,*) 'Correct data format of input matrix'
PAUSE 1001
1001 CONTINUE
STOP
C
C
999 WRITE(9,*) 'Sample matrix successfully prepared for
* split-weighted averaging.'
WRITE(9,*) 'Continue now with program
* Stats61.out (batch) or
* Stats62.out (batch)'
C
PAUSE 1000
1000 CONTINUE
STOP
END

6.9. Listing for Stats62.f

PROGRAM STATS
C
C Version 6.2, Michael Knappertsbusch,
C Original version created on 14.5.2001
C Modified 20.6.2002 by changing format to output file (4 instead of 3 digits for
C #of samples in Means + 95%ci files).
C Modified 10.11.2003: Instead ow writing to the outputfile the latest values for
C chamber number and coiling, this version inserts zeros and n at the respective chambers.
C (Stats does not calculate statistics for chamber number and coiling direction !)
C Modified from version 5.2 to 6.2, 19.12.2003: Adaptation to program KeelWidth93.out.
C
C Version 6.1 and 6.2 are exactly the same, except that version 6.2 writes the results
C to the output file using tabulators between columns.
C
C Program to calculate simple univariate statistics from morphological
C data that originate from outline analysis of foraminifera in keel view.
C The input data are arranged in matrices, that are stored in the extended
C output files, that were created with program KeelWidth93.out.
C
C !!! Note: Batch processing mode !!!
C
C Input files:
C There are two types of input files:
C 1.) A file with variable FILE_LIST (max. 20 chars long) containg the names of all
C individual files RESKW_sample_species_nChamb_COIL to be treated. The individual
C RESKW_sample_species_nChamb_COIL were first composed with program KeelWidth93.out.
C Note, that RESKW_sample_species_nChamb_COIL files must eventually be composed
C for chamber number and coling direction per species prior to feeding them to the
C Stats program.
C
C 2.) One or several files with the generalized name RESKW_sample_species_nChamb_COIL
C (in variable INPUT). The names of these RESKW_sample_species_nChamb_COIL files are
C always 31 chars long. They contain the age, specimen- and species names, number of
C points in the originally traced outlines, number of chambers in the last whorl,
C coiling direction, and the morphological measurements (X,Y,R,Area,D10,D90,D9010,
C DMIN,DMAX,PHI1,PHI2,ARESLO,ARESUP,DEV).
C
C Content of input files RESKW_sample_species_nChamb_COIL:
C
C Header information (the first 6 lines from the top)
C Data body including (from left to wright):
C SPECIMEN: Specimen
C AGE: Sample age in million years
C TPTS: Number of traced points
C SPECIES: Species name
C NCHAMB: Number of chambers in last whorl
C COIL: Coiling direction (s or d)
C X: Spiral height in µm
C Y: Shell diameter in keel view, in µm
C R: Ratio of X/Y, = flatness of shell in keel view
C Ar: Keel view area, in mm2
C D10: Lower keelwidth, at 10% of the axial extension, expressed as a fraction of X
C D90: Upper keelwidth, at 90% of the axial extension, expressed as a fraction of X
C D9010: Asolute difference of (D90-D10), expressed as fraction of the width of the shell [ABS(D90-D10)/DX]
C DMIN: DMIN=MIN(D10,D90)/DX
C DMAX: DMAX=MAX(D10,D90)/DX
C PHI1: Upper keel angle, in degrees
C PHI2: Lower keel angle, in degrees
C ARESLO: Lower keel residual area, in mm2
C ARESUP: Upper keel residual area, in mm2
C DEV: Deviation of total keel residual area from the value of AREA, in per cent
C
C
C
C Output is a single file called "Means + 95%ci" containing the age and names of the sample,
C the species names, chamber #, coiling, and the means and ±95% confidence intervals
C about the means for each of the numeric morphological characters.
C The extrema, variances and standard deviations, and upper and lower 95% confidence level
C are calculated but not written to the output file (can later be modified if needed).
C (In the batch version of this program this output file can be needed as to generate
C scatter plots for displaying the statistics).
C
C
C Explanation of other Variables:
C
C NMAX = Number of records in input file
C i = Number of numerical morphological variables
C X(i) = Values X1 through X14
C MINX(i),MAXX(i) = Minumum and maximum values of X1 through X14
C MEANX(i) = Means of X(i)
C VARX(i) = Variances of X(i)
C STDX(i) = Standard deviations of X(i)
C CIX(i) = ±95% confidence interval
C LCIX(i),UCIX(i) =Lower and upper 95% univariate
C confidence intervals for X(i) from the means
C SUMX(i) = Partial sums of X(i)
C
C
C Assignment of variables to morphological variables:
C X(1) = X
C X(2) = Y
C X(3) = R
C X(4) = Ar
C X(5) = D10
C X(6) = D90
C X(7) = D9010
C X(8) = DMIN
C X(9) = DMAX
C X(10) = PHI1
C X(11) = PHI2
C X(12) = ARESLO
C X(13) = ARESUP
C X(14) = DEV
C
C Name conventions (character variables):
C SPECIMEN = Specimen name (17 chars long)
C SPECIES = Species name (Code, 7 chars long)
C
C
C
C
INTEGER NMAX,TPTS
REAL X(1:14),MINX(1:14),MAXX(1:14)
REAL MEANX(1:14),VARX(1:14),STDX(1:14),CIX(1:14)
REAL LCIX(1:14),UCIX(1:14),SUMX(1:14)
REAL AGE,ROOT
CHARACTER*1 HEADER1,TAB
CHARACTER*1 COIL
CHARACTER*2 NCHAMB
CHARACTER*17 SPECIMEN
CHARACTER*7 SPECIES
CHARACTER*20 FILE_LIST
CHARACTER*31 INPUT
C
TAB=CHAR(9) ! Tabulator on a iMac (machine dependent value)
C
C
C
WRITE(9,*) '***************************************'
WRITE(9,*) '* *'
WRITE(9,*) '* Program Stats6.2 *'
WRITE(9,*) '* Analyzes keel view outline data from a matrix *'
WRITE(9,*) '* (batch processing mode) *'
WRITE(9,*) '* *'
WRITE(9,*) '* 19.12.2003 by Michael Knappertsbusch *'
WRITE(9,*) '* Only use in combination with KeelWidth93.out *'
WRITE(9,*) '* *'
WRITE(9,*) '****************************************'
C
WRITE(9,*) ' '
WRITE(9,*) ' '
WRITE(9,*) ' . . .REMEMBER. . .'
WRITE(9,*) 'first compose input files for #Ch and Coil'
WRITE(9,*) ' before running Stats6 on data !'
WRITE(9,*) ' '
C
C
C
C First, prepare the OUTPUT file:
OPEN(16,FILE='Means + 95%ci',STATUS='NEW')
C
C Write header line to output file:
C
WRITE(16,*) '*'
WRITE(16,89) TAB,TAB,TAB,TAB,TAB,TAB
1,TAB,TAB,TAB,TAB,TAB,TAB
2,TAB,TAB,TAB,TAB,TAB,TAB
3,TAB,TAB,TAB,TAB,TAB,TAB
4,TAB,TAB,TAB,TAB,TAB,TAB
5,TAB,TAB,TAB,TAB
C
89 FORMAT('Specimen',A1,'Age',A1
*,'#pts',A1,'Species',A1,'#Ch',A1
*,'Coil',A1,'#samples',A1,'X in µm'
*,A1,'±95ciX',A1,'Y in µm',A1,'±95ciY'
*,A1,'R=X/Y',A1,'±95%ciR',A1,'Ar in mm2'
*,A1,'±95%ciAR',A1,'D10%',A1,'±95%ciD10'
*,A1,'D90%',A1,'±95%ciD90',A1,'D9010%'
*,A1,'±95ciD9010',A1,'Dmin%',A1
*,'±95%ciDmin',A1,'Dmax%',A1
*,'±95%ciDmax',A1,'Phi1°',A1
*,'±95%ciPhi1',A1,'Phi2°',A1,'±95%ciPhi2'
*,A1,'ARESLO in mm2',A1,'±95%ciAreslo',A1
*,'ARESUP in mm2',A1,'±95%ciAresup',A1
*,'DEV in %',A1,'±95%ciDEV')
C
C
C
C
C Now, get list with input file names, open it, and read the name of an input file:
C
WRITE(9,*) 'Enter the name of the list
1 containg the file names (max 20 chars)'
READ(9,'(A20)') FILE_LIST
OPEN(17,FILE=FILE_LIST,STATUS='OLD')
C
C Main loop with samples starts here:
C
90 READ(17,'(A31)',END=2001) INPUT ! Read new file (new species) from FILE_LIST
WRITE(9,3) INPUT
3 FORMAT('. . .file ',A31,'. . .')
C
C Open INPUT and determine the number of specimens (record lines) in it:
C
OPEN(15,FILE=INPUT,STATUS='OLD') ! Input file
C
C
C Determine number of records in the input file.
C First, read header block again:
C
DO 4 I=1,6
READ(15,1) HEADER1
4 CONTINUE
C
NMAX=0
5 READ(15,8,END=1000) SPECIMEN
8 FORMAT(A17)
NMAX=NMAX+1
GOTO 5
1000 REWIND(15)
C
ROOT=SQRT(REAL(NMAX))
C
C
C Now proceed with statistical analysis:
C
C Read first 6 lines of header block:
DO 6 I=1,6
READ(15,1) HEADER1
6 CONTINUE
1 FORMAT(A1)
C
C
C
C Determine extrema:
C
C Read first data line:
READ(15,2) SPECIMEN,AGE,TPTS,SPECIES,
1NCHAMB,COIL,(X(I),I=1,14)
2 FORMAT(A17,2X,F5.2,2X,I4,2X,A7
1,1X,A2,2X,A1,2X,F6.1,2X,F6.1,2X
2,F6.4,2X,F7.4,2X,F7.2,2X,F7.2,2X
3,F7.2,2X,F7.2,2X,F7.2,2X,F7.3,2X
4,F7.3,2X,F7.4,7X,F7.4,6X,F7.3)
C
C
C
C Special case: Statistics if there is only one data record available (NMAX=1):
C
IF (NMAX.EQ.1) THEN
DO 7 I=1,14
MINX(I)=X(I)
MAXX(I)=X(I)
MEANX(I)=X(I)
VARX(I)=0.0
STDX(I)=0.0
CIX(I)=0.0
LCIX(I)=X(I)
UCIX(I)=X(I)
7 CONTINUE
GOTO 2000 ! Write results to output file and stop program.
END IF
C
C
C Now, proceed with calculation of extrema:
C
DO 10 I=1,14
MINX(I)=X(I)
MAXX(I)=X(I)
10 CONTINUE
C
C Read next data line:
15 READ(15,2,END=999) SPECIMEN,AGE,TPTS,
1 SPECIES,NCHAMB,COIL,(X(I),I=1,14)
C
DO 20 I=1,14
IF (X(I).LE.MINX(I)) THEN
MINX(I)=X(I)
END IF
C
IF (X(I).GT.MAXX(I)) THEN
MAXX(I)=X(I)
END IF
20 CONTINUE
GOTO 15
C
999 REWIND(15)
C
C
C
C
C
C Calculate means:
C
C Read header block again:
DO 25 I=1,6
READ(15,1) HEADER1
25 CONTINUE
C
C Initialize SUMX(i):
C
DO 30 I=1,14
SUMX(I)=0.0
30 CONTINUE
C
35 READ(15,2,END=888) SPECIMEN,AGE,TPTS,
1 SPECIES,NCHAMB,COIL,(X(I),I=1,14)
C
C Calculate partial sums:
C
DO 40 I=1,14
SUMX(I)=SUMX(I)+X(I)
40 CONTINUE
GOTO 35
C
888 DO 50 I=1,14
MEANX(I)=SUMX(I)/NMAX
50 CONTINUE
C
REWIND(15)
C
C
C
C
C Determine the variances, standard deviations and the 95% confidence intervals:
C
C Again, read header block:
DO 51 I=1,6
READ(15,1) HEADER1
51 CONTINUE
C
C Re-initialize SUMX(i):
C
DO 60 I=1,14
SUMX(I)=0.0
60 CONTINUE
C
C Calculate sums of squares of differences from the means:
C
65 READ(15,2,END=777) SPECIMEN,AGE,TPTS,
1 SPECIES,NCHAMB,COIL,(X(I),I=1,14)
C
DO 70 I=1,14
SUMX(I)=SUMX(I)+(X(I)-MEANX(I))**2
70 CONTINUE
GOTO 65
C
777 DO 80 I=1,14
VARX(I)=SUMX(I)/(NMAX-1)
STDX(I)=SQRT(VARX(I))
CIX(I)=1.96*STDX(I)/ROOT
LCIX(I)=MEANX(I)-CIX(I)
UCIX(I)=MEANX(I)+CIX(I)
80 CONTINUE
C
C
C
C Output of results to file "Means + 95%ci":
C
2000 WRITE(16,85) SPECIMEN(1:11),TAB,AGE,TAB,TPTS
*,TAB,SPECIES,TAB,' n',TAB,' C ',TAB,NMAX
*,TAB,MEANX(1),TAB,CIX(1),TAB,MEANX(2)
*,TAB,CIX(2),TAB,MEANX(3),TAB,CIX(3)
*,TAB,MEANX(4),TAB,CIX(4),TAB,MEANX(5)
*,TAB,CIX(5),TAB,MEANX(6),TAB,CIX(6)
*,TAB,MEANX(7),TAB,CIX(7),TAB,MEANX(8)
*,TAB,CIX(8),TAB,MEANX(9),TAB,CIX(9)
*,TAB,MEANX(10),TAB,CIX(10),TAB,MEANX(11)
*,TAB,CIX(11),TAB,MEANX(12),TAB,CIX(12)
*,TAB,MEANX(13),TAB,CIX(13),TAB,MEANX(14)
*,TAB,CIX(14)
C
85 FORMAT(A11,A1,F5.2,A1,I4,A1,A7,A1
1,A3,A1,A4,A1,I4,A1,F6.1,A1,F6.1
2,A1,F6.1,A1,F6.1,A1,F6.4,A1
3,F6.4,A1,F7.4,A1,F7.4,A1,F7.2,A1
4,F7.2,A1,F7.2,A1,F7.2,A1,F7.2,A1
5,F7.2,A1,F7.2,A1,F7.2,A1,F7.2,A1
6,F7.2,A1,F7.3,A1,F7.3,A1,F7.3,A1
7,F7.3,A1,F7.4,A1,F7.4,A1
8,F7.4,A1,F7.4,A1,F7.3,A1
9,F7.3)
C
CLOSE(15)
GOTO 90 ! Read new file from FILE_LIST
C
2001 PAUSE 100
100 CONTINUE
STOP
END

6.10. Listing for program Homolo31.f

PROGRAM HOMOL3
C
C Version 3.0, October 22, 1997 by M. Knappertsbusch
C
C !!!Batch processing!!!
C
C Modifications:
C 19.12.2000: Addition of an error report for the case that files are missing,
C and in the case that the specimen outline has not the same number of points
C than the reference outline.
C
C Correlation and homologization of the outlines of two digitized light objects.
C Homologization is performed by rotation in anti- and clockwise direction
C until the outlines match best (highest linear correlation coefficient).
C
C Input files are:
C
C - A file called Ref_char_nchamber_coil (21 characters long, interpreted by
C character variable REFERENCE containing the cartesian X and Y coordinates of
C a reference outline. "char" in the string is either Keel or Umbil for
C keel or umbilical view, respectively; n is the number of chambers;
C and coil gives the coining direction (sin or dex).
C
C - File with name FILE_LIST containing a list of the names of all files INPUT with
C the individual specimen outlines. FILE_LIST is a character variable of length 20.
C
C - All files, that contain the specimen outlines to be rotated and homologized with
C respect to the reference outline. The names of these files are read and interpreted
C through the character variable INPUT. The expression of INPUT is used to
C compose the outputfile with the x,y coordinates of the rotated and
C homologized specimen (via character variable OUTPUT) and the file with the
C correlation coefficients (via character variable CORRELATION).
C
C Outputfiles are:
C - Files ending with suffix _HOM: These are the cartesian x,y coordinates of
C the rotated and homologized specimens in position of the maximum correlation
C with respect to the reference specimen.
C The names of these files (interpreted through OUTPUT) are composed from the
C corresponding INPUT file.
C
C - Files ending with suffix _CORR: These files contain the linear correlation
C coefficients between rotated outlines of the specimen and the reference specimen
C as a function of the rotation angle. Given are the direction of the last rotational
C operation, the total angle of rotation (in degrees) and the correlation
C coefficients of the rotated speciem with respect to the reference outline for every
C position. The names (interpreted through CORRELATION) of these files are also
C composed from the corresponding INPUT file.
C
C
C Note: All file handling for morphometric analyses follows a name convention:
C
C filename = Name of sample, 15 characters long.
C filename_T = Input file for raw curves with cartesian coordinates, 17 characters long.
C filename_POL = Polar coordinates of filename, 19 characters long.
C filename_INT = Interpolated (cartesian) coordinates of filename, 19 characters long.
C filename_NORM = Normalized curve, 20 characters long.
C filename_HARM = Harmonic coefficients etc. of filename, 20 characters long.
C filename_AVG = Averaged harmonic coefficients over the entire sample 15 characters long.
C filename_REC = Reconstructed X,Y coordinates from fourier coefficients 15 characters long.
C coefficients
C filename_HOM = Cartesian coordinates of homologized specimen 19 characters long.
C filename_CORR = Linear correlation coefficient between refernce 20 characters long.
C specimen and the rotated specimen as a function of
C the rotation angle.
C
C
C Name convention for the reference outline (21 characters long):
C Ref_char_nChamber_coil char=Keel or Umbi (keel or umbilical view, respectively)
C n = chamber number in last whorl, as an integer number
C coil=coiling direction, sin or dex
C (for sinistral or dextral, respectively).
C
C
C
C
INTEGER I,M,MAX,NMAX,NSPEC,DIR
LOGICAL FILEX
DOUBLE PRECISION PI
PARAMETER(MAX=250, PI=3.14159265)
C
C*** Attention: When changing MAX adapt also MAX in subroutines MEANS and CORR !!!!
C
DOUBLE PRECISION RREF(1:MAX),RSPEC(1:MAX),R(1:MAX),PHI(1:MAX)
DOUBLE PRECISION SPEC(1:MAX),RAD(1:MAX),ALPHA(1:MAX)
DOUBLE PRECISION XREF,YREF,XSPEC,YSPEC,RHO,THETA
DOUBLE PRECISION MREF,MSPEC,MEAN,XHOM,YHOM
DOUBLE PRECISION DTHETA,COEF,MCOEF,ROT
CHARACTER*11 HEADER
CHARACTER*19 OUTPUT
CHARACTER*20 INPUT,CORRELATION,FILE_LIST
CHARACTER*21 REFERENCE

C
C Explanation of variables:
C
C Vector RREF contains the radii of the reference outline.
C Vector RSPEC contains the radii of the specimen outline.
C Vector SPEC contains the unrotated radii of the specimen outline.
C Vector PHI contains the angular arguments of the unrotated specimen outline (in radians).
C Vector R contains the locally rotated radii of the specimen.
C Vector RAD contains the radii of the specimen outline in homologous position to the reference
C outline.
C Vector ALPHA contains the angular arguments of the specimen outline in homologous position
C (in radians).
C Variables XREF and YREF are used to read the cartesian coordinates of the reference outline.
C Variables XSPEC and YSPEC are used to read the cartesian coordinates of the specimen outline.
C Variables RHO and THETA (in radians) are locally used to convert cartesian into polar coordinates.
C Variable DTHETA is the angular increment when the specimen is rotated (in radians).
C Variable COEF is the local correlation coefficient between the radii of the rotated specimen against
C the reference outline.
C Variable MCOEF is the maximum correlation coefficient reach during the rotations (= homologous
C position).
C Variable ROT is the rotation angle (in degrees !) between the reference outline and the specimen
C outline at homologous position.
C
WRITE(9,*) 'Homologization of specimens to
* a reference specimen'
WRITE(9,*) '. . .Enter name of file with
* reference outline (21 chars)'
READ(9,4) REFERENCE
4 FORMAT(A21)
C
C Open file with reference specimen and keep it open:
C
OPEN(15,FILE=REFERENCE,STATUS='OLD')
C
C
C Recording errors during program run:
C
OPEN(25,FILE='Errors_during_Homolo',STATUS='NEW')
C
C
WRITE(9,*) '. . .Names of interpolated and normalized
* outline data needed. . .'
WRITE(9,*) ' Enter list containing files (20 chars)'
READ(9,1) FILE_LIST
1 FORMAT(A20)
C
C Now opening that file list and reading a name of the file containing the x,y coordinates of
C an interpolated and normalized specimen outline:
C
OPEN(14,FILE=FILE_LIST,STATUS='OLD')
C
C Determine the sample name from INPUT, and write the sample name into the
C header of file Errors_during_Homolo:
C
READ(14,2) INPUT
HEADER=INPUT(1:11)
WRITE(25,7) HEADER
7 FORMAT('Sample name: ',A11)
REWIND(14)
C
C
C*** Main loop starts at statement #5:
C
5 READ(14,2,END=999) INPUT
2 FORMAT(A20)
WRITE(9,3) 'Specimen ',INPUT
3 FORMAT(A9,A20)
C
C
C*** Check first, whether file with name INPUT exists or not. If it does not exist:
C print a message to the error log 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
8 FORMAT('File ',A20,' is missing')
WRITE(9,*) '. . .Program continues. . .'
GOTO 5
END IF
C
C
C Next, compose names for output file with correlation coefficients
C and file with X,Y coordinates of rotated and homologized specimen from
C the name of the input file INPUT:
C
CORRELATION=INPUT(1:15)//'_CORR'
OUTPUT=INPUT(1:15)//'_HOM'
C
C
C
C Now, reading reference outline:
C NMAX is the number of points in the reference outline.
C

I=1
10 READ(15,*,END=20) XREF,YREF
C
C Convert into polar coordinates:
C
CALL POLAR(XREF,YREF,RHO,THETA)
RREF(I)=RHO
I=I+1
GOTO 10
C
C
20 REWIND(15)
C
C Determine number of points in reference file:
C
NMAX=I-1
C
IF (NMAX.GT.MAX) THEN
WRITE(9,*) 'Adapt dimensions in program'
PAUSE 21
21 CONTINUE
STOP
END IF
C
C Calculate angular increments:
C
DTHETA=2*PI/NMAX

C
C Preparation for the calculation of correlation coefficients between reference outline
C and specimen outline: Calculation of the mean MREF of the components of the
C vector RREF(I).
C
CALL MEANS(RREF,NMAX,MEAN)
MREF=MEAN
C
C
C Now, reading outline of specimen to be measured:
C
OPEN(16,FILE=INPUT,STATUS='OLD')
I=1
30 READ(16,*,END=40) XSPEC,YSPEC
C
C Convert into polar coordinates:
C
CALL POLAR(XSPEC,YSPEC,RHO,THETA)
RSPEC(I)=RHO
PHI(I)=THETA
I=I+1
GOTO 30
C
C
C*** Check whether specimen outline has the same number of points as the reference outline:
C
40 NSPEC=I-1
IF (NSPEC.NE.NMAX) THEN
WRITE(9,*) 'Reference and specimen outlines do not have
* the same number of points !!!'
WRITE(9,*) '(Reference outline has ',NMAX,' points)'
WRITE(9,*) '(Specimen outline has ',NSPEC,' points)'
WRITE(9,*) '. . .Program reads next sample. . .'
WRITE(25,41) INPUT,NSPEC,NMAX
41 FORMAT('Error: ',A20,': has ',I4,' points,
* Ref outline has ',I4,' points')
GOTO 5
END IF
CLOSE(16)
C
C
C Memorize radii of specimen outline in separate vector SPEC(I):
C (The unrotated specimen or SPEC will be used later for the rotation
C in clockwise direction).
C
DO 45, I=1,NMAX
SPEC(I)=RSPEC(I)
45 CONTINUE
C
C
C Again preparation for the calculation of correlation coefficients between reference outline
C and specimen outline: Calculation of the mean MSPEC of the components of the
C vectors RSPEC(I).
C
CALL MEANS(RSPEC,NMAX,MEAN)
MSPEC=MEAN
C
C
C Now, calculate linear correlation coefficients of the values of RREF and
C the unrotated RSPEC, and write results to file "HOMOLOGIZED":
C
OPEN(17,FILE=CORRELATION,STATUS='NEW')
C
C Writing header:
WRITE(17,*) 'Last direction, total direction
* (in degrees), correlation coefficient'
WRITE(17,*) ' '
WRITE(17,*) '1. Rotation in anti-clockwise direction):'
C
C The integer variable M counts how many times the specimen outline has been rotated.
C The integer variable DIR indicates the direction of the latest completed rotation; a positive
C value of DIR means a rotation in anti-clockwise direction, and a negative value means a
C rotation in a clockwise direction.
C
M=0
DIR=0
CALL CORR(RREF,RSPEC,NMAX,MREF,MSPEC,COEF)
C
C Memorize first correlation coefficient:
C
MCOEF=COEF
C
C Write result to file:
C
WRITE(17,*) DIR,',',M,',',MCOEF
C
C Now, rotation of specimen outline by increment DTHETA and calculation of
C new correlation coefficient. Rotated outline is stored in vector R.
C First rotation in anti-clockwise direction (= +DIR), if necessary, subsequent
C rotation in clockwise direction (= -DIR).
C
C Rotation in anti-clockwise direction:
C First increase M by one to indicate the next rotational operation.
C
70 M=M+1
DIR=1
DO 50, I=1,NMAX-1
R(I)=RSPEC(I+1)
50 CONTINUE
C
C Endpunkt rotieren:
C
R(NMAX)=RSPEC(1)
C
C Calculate new mean of components of the rotated specimen outline:
C
CALL MEANS(R,NMAX,MEAN)
MSPEC=MEAN
C
C Calculation of new correlation coefficient:
C
CALL CORR(RREF,R,NMAX,MREF,MSPEC,COEF)
C
C Find position (=ROT) with maximum correlation coefficient (=MCOEF).
C Store polar coordinates in this position in vectores RAD and ALPHA.
C
IF (COEF.GE.MCOEF) THEN
MCOEF=COEF
ROT=M*DTHETA*180/PI
DO 55, I=1,NMAX
RAD(I)=SPEC(I)
ALPHA(I)=PHI(I)+M*DTHETA
55 CONTINUE
END IF
C
WRITE(17,*) DIR,',',M*DTHETA*180/PI,',',COEF
C
C Stop at +30 degrees:
C
IF (M*DTHETA.GE.PI/6) THEN
WRITE(9,*) '. . .Rotation has reached +30 degrees. . .'
GOTO 600
END IF
C
C Next rotational cycle:
C
C First initialize next rotational cycle
C by replacing previous RSPEC(I) by new R(I):
C
DO 65, I=1,NMAX
RSPEC(I)=R(I)
65 CONTINUE
C
GOTO 70
C
C
C
C
C Rotation of specimen outline in clockwise direction, and recalculation of correlation
C coefficient.
600 WRITE(17,*) ' '
WRITE(17,*) '2. Rotation in clockwise direction:'
WRITE(17,*) ' '
C
C First, reset the specimen in the unrotated position.
C
DO 80, I=1,NMAX
RSPEC(I)=SPEC(I)
80 CONTINUE
C
C
C Now, calculate again linear correlation coefficients of the values of RREF and
C the unrotated RSPEC, and write results to the file with the correlation
C coefficients (CORRELATIONS):
C
M=0
DIR=0
CALL CORR(RREF,RSPEC,NMAX,MREF,MSPEC,COEF)
C
C Check again for position with maximum correlation coefficient, and
C store coordinates in vectors RAD and ALPHA:
C
IF (COEF.GE.MCOEF) THEN
MCOEF=COEF
ROT=M*DTHETA*180/PI
DO 85, I=1,NMAX
RAD(I)=SPEC(I)
ALPHA(I)=PHI(I)
85 CONTINUE
END IF
C
C Write result to file:
C
WRITE(17,*) DIR,',',M,',',COEF
C
C Now, rotation in clockwise direction:
C
90 M=M+1
DIR=-1
R(1)=RSPEC(NMAX)
DO 100, I=2,NMAX
R(I)=RSPEC(I-1)
100 CONTINUE
C
C Calculation of new correlation coefficient:
C
CALL CORR(RREF,R,NMAX,MREF,MSPEC,COEF)
C
C Again, find position with maximum correlation coefficient, and
C store rotated specimen in best position in vectors RAD and ALPHA:
C
IF (COEF.GE.MCOEF) THEN
MCOEF=COEF
ROT=-M*DTHETA*180/PI
DO 102, I=1,NMAX
RAD(I)=SPEC(I)
ALPHA(I)=PHI(I)-M*DTHETA
102 CONTINUE
END IF
C
WRITE(17,*) DIR,',',-M*DTHETA*180/PI,',',COEF
C
C Stop at -30 degrees, and output of results:
C
IF (M*DTHETA.GE.PI/6) THEN
WRITE(9,*) '. . .Rotation has reached -30 degrees. . .'
WRITE(9,105) 'The maximum correlation (',MCOEF,') is
* at an angle of ',ROT, ' degrees.'
105 FORMAT(A25,F5.3,A20,F5.1,A9)
C
C Write cartesian coordinates (XHOM,YHOM) of homologized specimen
C to file OUTPUT:
C
OPEN(18,FILE=OUTPUT,STATUS='NEW')
C
C Re-ronvert into cartesian coordinates:
C
DO 107, I=1,NMAX
XHOM=RAD(I)*COS(ALPHA(I))
YHOM=RAD(I)*SIN(ALPHA(I))
WRITE(18,*) XHOM,',',YHOM
107 CONTINUE
C
C Read next specimen from file list:
C
GOTO 5
C
999 WRITE(9,*) 'All specimens treated'
CLOSE(25)
C
WRITE(9,*) 'Enter a key'
PAUSE789
789 CONTINUE
STOP
END IF
C
C Next rotational cycle:
C
C First initialize next rotational cycle
C by replacing previous RSPEC(I) by new R(I):
C
DO 110, I=1,NMAX
RSPEC(I)=R(I)
110 CONTINUE
C
GOTO 90
C
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=2*PI/3
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

SUBROUTINE MEANS(VEC,KMAX,MEA)
C
C Subroutine to calculate the mean (= MEA) from the components of the one-dimensional
C vector VEC. KMAX indicates the actual dimension of VEC, with KMAX<=MAX.
C Input into the subroutine is VEC and KMAX.
C The subroutine returns the value of MEA to the calling unit.
C
INTEGER MAX,KMAX,I
PARAMETER(MAX=250)
DOUBLE PRECISION VEC(1:MAX)
DOUBLE PRECISION MEA,SUM
C
SUM=0.0
DO 10, I=1,KMAX
SUM=SUM+VEC(I)
10 CONTINUE
C
MEA=SUM/KMAX
C
RETURN
END

SUBROUTINE CORR(REF,SPEC,KMAX,MEREF,MESPEC,CORREL)
C
C Subroutine to calculate the linear correlation coefficient (= CORREL) of the components
C of REF and SPEC. Input into the subroutine are the two one-dimensional vectors
C REF and SPEC containing KMAX components each, and KMAX, which must be smaller or
C equal to MAX, the maximum dimension of REF and SPEC. MESPEC and MEREF are also
C input values and are the means of the components of SPEC and REF, respectively.
C Output is the single variable CORREL, which is returned to the calling unit.
C
INTEGER KMAX,MAX,I
PARAMETER(MAX=250)
DOUBLE PRECISION REF(1:MAX),SPEC(1:MAX)
DOUBLE PRECISION MESPEC,MEREF,CORREL
DOUBLE PRECISION SUM1,SUM2,SUM3
C
C Initialization of individual sums:
C
SUM1=0.0
SUM2=0.0
SUM3=0.0
C
C Calculate individual sums:
C
DO 10, I=1,KMAX
SUM1=SUM1+(REF(I)-MEREF)*(SPEC(I)-MESPEC)
SUM2=SUM2+(REF(I)-MEREF)**2
SUM3=SUM3+(SPEC(I)-MESPEC)**2
10 CONTINUE
C
C Calculation of the correlation coefficient:
C
CORREL=SUM1/SQRT(SUM2*SUM3)
C
RETURN
END

6.11. Listing for Program Averou3.f

PROGRAM AVEROU
C
C Version 3.0, June 21, 2002 (with option for split-weighted averages).
C Version 2.0, October 22, 1997, by M. Knappertsbusch.
C
C !!! Processing in batch mode. !!!
C
C Given a series of outlines from a digitized light object this program calculates an average
C outline from the cartesian X and Y coordinates.
C Attention: All outlines must be in homologous position !!!
C
C There are 2 options:
C - Without split-weighted averages (=option 1)
C - With split-weighted averages (=option 2)
C
C Split-weighted averaging is necessary if in a single sample foraminiferal
C specimens were unevenly distributed in different size classes.
C
C Example: Only 2 specimens of G. menardii "A" occur in entire split of size
C class 500-1000µm, but 50 specimens of the same species occur in 1/32 split of
C size class 100-500µm.
C
C In this example, the statistical influence of the larger size fraction is small,
C and the averaged outline is mainly made by the smaller specimens. Option 2
C accounts for this effect, so that in the above example each of the two specimens
C from the larger size fraction contribute little (1/1 split), but each specimen of
C the smaller size fraction contributes 32 times as much to the average outline.
C
C
C Input files are:
C - One file (=FILE_LIST) containing a list of the file names (=INPUT) where the cartesian
C X and Y coordinates of the individual outlines are written, from which the average outline
C is calculated. FILE_LIST is a character variable of 20 characters lengths.
C
C Note:
C In option 1, FILE_LIST contains only the names, i.e. Format(A19)
C In option 2, FILE_LIST has a different format, with the first two digits are
C reserved for the plit number SPLIT (i.e. the number of times a file must be repeated
C in FILE_LIST_spw), then the size fraction is written in 10 characters, and then the
C names (19 chars) of the individual outlines are written. The format for option 2 is
C FORMAT(I2,2X,A10,2X,A19).
C
C - All files, that contain the X,Y (cartesian) coordinates of the individual
C outlines. The names of these files are interpreted through the character
C variable INPUT. The expression of INPUT is used to compose the name of the
C output file (suffix _AVOUTLIN).
C
C Output files are:
C - One single file with the suffix _AVOULIN containing the cartesian X,Y coordinates of the
C desired mean outline.
C
C
C Note: All file handling for morphometric analyses follows a name convention:
C
C filename = Name of sample, 15 characters long.
C filename_T = Input file for raw curves with cartesian coordinates, 17 characters long.
C filename_POL = Polar coordinates of filename, 19 characters long.
C filename_INT = Interpolated (cartesian) coordinates of filename, 19 characters long.
C filename_NORM = Normalized curve, 20 characters long.
C filename_HARM = Harmonic coefficients etc. of filename, 20 characters long.
C filename_AVG = Averaged harmonic coefficients over the entire sample 15 characters long.
C filename_REC = Reconstructed X,Y coordinates from fourier coefficients 15 characters long.
C coefficients
C filename_HOM = Cartesian coordinates of homologized specimen 19 characters long.
C filename_CORR = Linear correlation coefficient between refernce 20 characters long.
C specimen and the rotated specimen as a function of
C the rotation angle.
C filename_AVOULIN = Mean outline obtained through 1.) rotation and homologiza-20 characters long.
C tion of individual outlines, 2.) averaging outlines with
C program Averou.f.
C
C
C Name convention for the reference outline (21 characters long):
C Ref_char_nChamber_coil char=Keel or Umbi (keel or umbilical view, respectively)
C n = chamber number in last whorl, as an integer number
C coil=coiling direction, sin or dex
C (for sinistral or dextral, respectively).
C
C
C
INTEGER I,MAX,NMAX,N,SPLIT
PARAMETER(MAX=250)
DOUBLE PRECISION SRAD(1:MAX),ANGLE(1:MAX)
DOUBLE PRECISION MRAD(1:MAX)
DOUBLE PRECISION NRAD(1:MAX)
DOUBLE PRECISION X,Y,R,THETA,XMEAN,YMEAN
CHARACTER*19 INPUT
CHARACTER*20 FILE_LIST,OUTPUT
CHARACTER*1 ANSWER
CHARACTER*10 SFRACT
C
C Explanation of variables:
C MAX determines the maximum dimension of the vectors
C SRAD,ANGLE,MRAD and NRAD.
C N is the number of outlines to be averaged.
C NMAX is the number of points in each individual outline
C (note: all outlines must have the same number of points).
C
C X,Y and R,THETA are the cartesian and polar coordinates of a single points of
C an outline, respectively.
C SRAD and SANGLE are the partial sums of points of subsequently read outlines. They are
C one-dimensional vectors and their components are in polar coordinates.
C NRAD are the radii of the next outline (polar coordinates, at similar angular arguments).
C MRAD contains the radii of the average outline.
C
WRITE(9,*) 'Calculation of average outline'
C
WRITE(9,*) 'Are all outlines in a homologous position ?'
WRITE(9,*) ' y/n'
READ(9,2) ANSWER
2 FORMAT(A1)
IF (ANSWER.EQ.'n') THEN
STOP
END IF
C
C
C Choose option (1=no split-weighted averaging of outlines,
C 2=with split-weighted averaging of outlines):
C
WRITE(9,*) 'Simple averages (=1) or
* split-weighted averages (=2) ?'
C
WRITE(9,*) ' '
WRITE(9,*) 'Note different input formats in FILE_LIST:'
WRITE(9,*) 'Option 1: FORMAT(A19)'
WRITE(9,*) 'Option 2: FORMAT(I2,2X,A10,2X,A19)'
WRITE(9,*) ' '
C
READ(9,2) ANSWER
C
C
C
C Get name of the file list with all filenames to be treated:
C
WRITE(9,*) ' . . .Names of outlines needed. . .'
WRITE(9,*) 'enter list containing all files (20 chars)'
READ(9,3) FILE_LIST
3 FORMAT(A20)
C
C
C In case of option 2: Read the split factor SPLIT from FILE_LIST, open a new file FILE_LIST_spw,
C and write the names of the INPUT files according to the factor SPLIT. Subsequently use
C FILE_LIST_spw as the new list with input files.
C
IF (ANSWER.EQ.'2') THEN
OPEN(12,FILE='FILE_LIST_spw',STATUS='NEW')
OPEN(14,FILE=FILE_LIST,STATUS='OLD')
8 READ(14,6,END=9) SPLIT,SFRACT,INPUT
6 FORMAT(I2,2X,A10,2X,A19)
DO 10, I=1,SPLIT
WRITE(12,7) INPUT
10 CONTINUE
7 FORMAT(A19)
GOTO 8
END IF
9 CLOSE(12)
CLOSE(14)
C
C
C Now, open the file list (Note: Individual file names in FILE_LIST must be not more
C than 19 characters long). Read first file name in it, open that file and read cartesian X,Y
C coordinates:
C
IF (ANSWER.EQ.'2') THEN
FILE_LIST='FILE_LIST_spw'
END IF
C
OPEN(15,FILE=FILE_LIST,STATUS='OLD')
READ(15,5) INPUT
write(9,5) INPUT
5 FORMAT(A19)
C
C Compose name of output file:
C
OUTPUT=INPUT(1:11)//'_AVOUTLIN'
C
C Initialization of counter for number of outlines:
C
N=1
C
OPEN(16,FILE=INPUT,STATUS='OLD')
I=0
20 I=I+1
READ(16,*,END=99) X,Y
C
C Calculate polar coordinates and store in vectors SRAD and SANGLE (= initialization for
C partial sum and associated angular argument).
C
CALL POLAR(X,Y,R,THETA)
SRAD(I)=R
ANGLE(I)=THETA
C
IF (I.GT.MAX) THEN
WRITE(9,*) 'Adapt dimensions in program !'
PAUSE 19
19 CONTINUE
STOP
END IF
C
GOTO 20
C
C Determine number of points in outline:
C
99 NMAX=I-1
CLOSE(16)
C
C Read next file name open file and read X,Y coordinates:
C
25 READ(15,5,END=9999) INPUT
C
C Increase counter N for number of outlines by one:
C
N=N+1
C
OPEN(17,FILE=INPUT,STATUS='OLD')
I=0
30 I=I+1
READ(17,*,END=999) X,Y
CALL POLAR(X,Y,R,THETA)
NRAD(I)=R
GOTO 30
C
999 CLOSE(17)
C
C Now, calculate partial sums for all radii (angular arguments remain the same !!!):
C
DO 40, I=1,NMAX
SRAD(I)=SRAD(I)+NRAD(I)
40 CONTINUE
C
C Read next filename, open it, read it and calculate next partial sums:
C
GOTO 25
C
C Calculate average outline, re-convert into cartesian coordinates, open
C output file and write results to output file:
C
9999 DO 50, I=1,NMAX
MRAD(I)=SRAD(I)/N
50 CONTINUE
C
OPEN(18,FILE=OUTPUT,STATUS='NEW')
DO 60, I=1,NMAX
XMEAN=MRAD(I)*COS(ANGLE(I))
YMEAN=MRAD(I)*SIN(ANGLE(I))
WRITE(18,*) XMEAN,',',YMEAN
60 CONTINUE
C
WRITE(9,*) '. . .One file ',OUTPUT,' written'
61 FORMAT(A14,A20,A8)
WRITE(9,*) '. . .Number of outlines averaged: ',N
WRITE(9,*) 'Number of points in averaged outline: ',NMAX
PAUSE 500
500 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=2*PI/3
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

6.12. Listing for program Harman.f

PROGRAM HARMAN
C
C Version 2, October 13, 1997, by Michael Knappertsbusch
C Minor corrections added July 3, 2005.
C
C !!! Processing in batch mode !!!
C
C Harmonic analysis of a closed plane curve consisting of
C N discrete points. Points are given as cartesian coordinates
C U and V. The program converts U and V into polar coordinates
C R and Theta and then goes through a polar fourier analysis.
C The number of harmonics (K) to be calculated is determined by KMAX.
C Output are the coefficients ALPHAK and BETAK of the Fourier
C series as a function of K, the corresponding coefficients and
C phase angle AK and PHIK of the polar form of the fourier series,
C and the power POWERK of the discrete spectrum of the curve as
C a function of K.
C
C Parameters are calculated after Davis, 1986, pages 248 ff.
C
C Input files are:
C - File with name FILE_LIST containing a list of all
C names of files (filename), that contain the cartesian
C coordinates of the interpolated and normalized outlines.
C FILE_LIST is a character variable and can bear any
C name of 20 characters length.
C
C - All files, that contain the X,Y coordinates of the
C outlines, that have previously been interpolated and
C normalized; i.e. all files ending with suffix _NORM).
C The names of these files are interpreted through the
C character variable INPUT. The expression of INPUT is
C used to compose the name of the corresponding output
C file (suffix _HARM), which contains the harmonics
C coefficients, phase angles and the power of the spectra.
C
C Output files are:
C - One file per outline, that contains the harmonics
C coefficients, phase angles and the power of the spectra.
C The names of the output files are also composed through
C the character variable INPUT by addition of the suffix
C _HARM to the name of filename.
C
C Note: All file handling for morphometric analyses follows a name convention:
C
C filename = Name of sample, 15 characters long.
C filename_T = Input file for raw curves with cartesian coordinates, 17 characters long.
C filename_POL = Polar coordinates of filename, 19 characters long.
C filename_INT = Interpolated (cartesian) coordinates of filename, 19 characters long.
C filename_NORM = Normalized curve, 20 characters long.
C filename_HARM = Harmonic coefficients etc. of filename, 20 characters long.
C
INTEGER K,KMAX,N
DOUBLE PRECISION PI
DOUBLE PRECISION SUMR,SUMA,SUMB,RJ,THETAJ
DOUBLE PRECISION U,V,ALPHAK,BETA0,BETAK,AK,PHIK
DOUBLE PRECISION POWERK
CHARACTER*20 FILE_LIST,OUTPUT
CHARACTER*20 INPUT
C
PARAMETER(PI=3.14159265)
C
WRITE(9,*) '********************************************'
WRITE(9,*) '* *'
WRITE(9,*) '* Harmonic decomposition of interpolated *'
WRITE(9,*) '* and normalized outlines *'
WRITE(9,*) '* *'
WRITE(9,*) '* Processing in batch mode *'
WRITE(9,*) '* *'
WRITE(9,*) '* Version 2.1, 3 July 2005 *'
WRITE(9,*) '* By Michael Knappertsbusch *'
WRITE(9,*) '* *'
WRITE(9,*) '********************************************'
WRITE(9,*) ' '
C WRITE(9,*) 'Enter the number of harmonics to calculate (KMAX)'
C READ(9,*) KMAX

WRITE(9,*) '. . .Names of interpolated
1and normalized outlines needed. . .'
WRITE(9,*) ' enter list containing names
2 of files (20 chars)'
READ(9,2) FILE_LIST
2 FORMAT(A20)
C
OPEN(20,FILE=FILE_LIST,STATUS='OLD')
4 READ(20,6,END=999) INPUT
6 FORMAT(A20)
WRITE(9,3) INPUT
3 FORMAT('. . .Reading new file ',A20,'. . .')
C
C Determination of actual filenames:
C
OUTPUT=INPUT(1:15)//'_HARM'
C
OPEN(18,FILE=INPUT,STATUS='OLD')
OPEN(19,FILE=OUTPUT,STATUS='NEW')
C
C
N=O
SUMR=0.0
C
C***** First calculation of zero harmonic (K=0):
C
K=0
C
5 READ(18,*,END=10) U,V
N=N+1
C
C Transform cartesian coordinates into polar coordinates RJ
C and THETAJ:
C
CALL POLAR(U,V,RJ,THETAJ)
C
SUMR=SUMR+RJ
GOTO 5
C
C Calculate last sum:
C
10 CALL POLAR(U,V,RJ,THETAJ)
SUMR=SUMR+RJ
N=N+1
C
C Determine necessary number of Kmax:
C
KMAX=INT(N-1)/2-1
C
ALPHAK=0.0
BETA0=SUMR/N
AK=BETA0
PHIK=0.0
POWERK=(AK**2)/2
C
WRITE(19,*) K,',',ALPHAK,',',BETA0,',',AK,',',PHIK,',',POWERK
C
REWIND(18)
C
C
C***** Now calculation of all other harmonics:
C
DO 15, K=1,KMAX
SUMA=0.0
SUMB=0.0
20 READ(18,*,END=25) U,V
CALL POLAR(U,V,RJ,THETAJ)
SUMA=SUMA+RJ*COS(K*THETAJ)
SUMB=SUMB+RJ*SIN(K*THETAJ)
GOTO 20
C
C Last sum:
C
25 CALL POLAR(U,V,RJ,THETAJ)
SUMA=SUMA+RJ*COS(K*THETAJ)
SUMB=SUMB+RJ*SIN(K*THETAJ)
C
ALPHAK=2*SUMA/N
BETAK=2*SUMB/N
AK=SQRT(ALPHAK**2+BETAK**2)
PHIK=ARCTAN(ALPHAK,BETAK)
POWERK=(AK**2)/2
C
WRITE(19,*) K,',',ALPHAK,',',BETAK,',',AK,',',PHIK,',',POWERK
C
REWIND (18)
15 CONTINUE
C
GOTO 4
C
999 WRITE(9,*) 'Analysis for all curves terminated'
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=2*PI/3
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***************************************************************************
DOUBLE PRECISION FUNCTION ARCTAN(US,VS)
C***************************************************************************
C
C Calculates the arcus tangent as a function of the X-coordinate US
C and the Y-coordinate VS. All angular values are in Radians.
C ARCTAN is returned in Radians.
C
DOUBLE PRECISION US,VS
DOUBLE PRECISION PI
C
PARAMETER(PI=3.14159265)
C
IF (US.GT.0.AND.VS.EQ.0) THEN
ARCTAN=0
ELSE IF (US.EQ.0.AND.VS.GT.0) THEN
ARCTAN=PI/2
ELSE IF (US.LT.0.AND.VS.EQ.0) THEN
ARCTAN=PI
ELSE IF (US.EQ.0.AND.VS.LT.0) THEN
ARCTAN=2*PI/3
END IF
C
C For all other cases:
C
IF (US.GT.0.AND.VS.GT.0) THEN
ARG=ABS(VS/US)
ARCTAN=ATAN(ARG)
ELSE IF (US.LT.0.AND.VS.GT.0) THEN
ARG=ABS(US/VS)
ARCTAN=ATAN(ARG)+PI/2
ELSE IF (US.LT.0.AND.VS.LT.0) THEN
ARG=ABS(VS/US)
ARCTAN=ATAN(ARG)+PI
ELSE IF (US.GT.0.AND.VS.LT.0) THEN
ARG=ABS(US/VS)
ARCTAN=ATAN(ARG)+3*PI/2
END IF
C
RETURN
END

6.13. Listing for program Recon2.f

PROGRAM RECON2
C
C Version 2, 22 October 1997, by Michael Knappertsbusch
C Method after Ehrlich and Weinberg
C
C Program to reconstruct a discretized curve from the coefficients
C of its discrete fourier approximation (K harmonics).
C
C Note: All angular values are in Radians.
C
C
C Input files are:
C - One file containing the U,V (cartesian) coordinates of the
C interpolated and normalized original curve. This filename is
C interpreted through the character variable INPUT1.
C - A second file containing the fourier coefficients for the K harmonic
C functions. This file is interpreted through the character variable INPUT2.
C
C Output files are:
C - One single file filename_REC, which is treated through the character
C variable OUTPUT. The value of output is composed through the value of
C INPUT2: It takes the first 11 characters of INPUT2 and adds the
C suffix _REC.
C
C
C Note: All file handling for morphometric analyses follows a name convention:
C
C filename = Name of sample, 15 characters long.
C filename_T = Input file for raw curves with cartesian coordinates, 17 characters long.
C filename_POL = Polar coordinates of filename, 19 characters long.
C filename_INT = Interpolated (cartesian) coordinates of filename, 19 characters long.
C filename_NORM = Normalized curve, 20 characters long.
C filename_HARM = Harmonic coefficients etc. of filename, 20 characters long.
C filename_AVG = Averaged harmonic coefficients over the entire sample 15 characters long.
C filename_REC = Reconstructed X,Y coordinates from fourier coefficients 15 characters long.
C coefficients
C
INTEGER K,KMAX
DOUBLE PRECISION U,V,R,THETA,SUM
DOUBLE PRECISION AK,BK,RK,PHIK,POWERK
DOUBLE PRECISION X,Y
CHARACTER*20 INPUT1,INPUT2
CHARACTER*15 OUTPUT
C
C
C
WRITE(9,*) '. . .Reconstruction of shape. . .'
WRITE(9,*) ' '
WRITE(9,*) 'Enter filename with normalized or interpolated data'
WRITE(9,*) ' (max. 20 characters'
C
C Reconstruction of shape can be done on either the normalized or interpolated data:
C Only the angle THETA will be used in the calculations.
C
READ(9,1) INPUT1
WRITE(9,*) 'Enter filename with the harmonic coefficients'
WRITE(9,*) ' (max 20 characters)'
READ(9,1) INPUT2
1 FORMAT(A20)
C
C Composition of the name of the output file from the name of
C the file containing the harmonic coefficients:
C
OUTPUT=INPUT2(1:11)//'_REC'
C
OPEN(15,FILE=INPUT1,STATUS='OLD')
OPEN(16,FILE=INPUT2,STATUS='OLD')
OPEN(17,FILE=OUTPUT,STATUS='NEW')
C
C Read the cartesian U and V coordinates U,V from
C file containing the interpolated values of the curve:
C
WRITE(9,2) '...Reconstructing from ',INPUT2
2 FORMAT(A23,A20)
WRITE(9,*) '...READING ',INPUT1,'...'
WRITE(9,3) '...Writing results to ',OUTPUT
3 FORMAT(A22,A15)
C
C
C
C
5 READ(15,*,END=999) U,V
C
C Convert U and V into polar coordinates R,THETA:
C
CALL POLAR(U,V,R,THETA)
C
C Now calculate local radii as the sum of harmonic terms at angle THETA by
C adding the periodic functions for each harmonic number:
C Read first R0:
C
READ(16,*) K,AK,BK,R0,PHIK,POWERK
SUM=R0
C
C Now adding all subsequent trigonometric terms:
C
20 READ(16,*,END=25) K,AK,BK,RK,PHIK,POWERK
SUM=SUM+AK*COS(K*THETA)+BK*SIN(K*THETA)
GOTO 20
C
C Last sum:
25 SUM=SUM+AK*COS(K*THETA)+BK*SIN(K*THETA)
C
C Reconvert into cartesian coordinates:
C
X=SUM*COS(THETA)
Y=SUM*SIN(THETA)
C
WRITE(17,*) X,',',Y
REWIND(16)
GOTO 5
C
999 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=2*PI/3
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

6.14. Listing for program Rename_Raw.f

program rename_raw
c
c Program to rename raw image files.
c By Michael Knappertsbusch, 11.8.2004.
c
c The program reads the contents of the files and writes them under a new name to disk.
c
c
c Explanations:
c Raw images have size of 640x480 pixels, 8 bit depth (256 grey levels).
c A raw image has no header. But it is useful (for programming purpose) to apply
c "word lengths" that correspond to the length of a Tiff header:
c By standard the first 768 bytes of a Tiff file belong to the Tiff header (according
c to Documentation of Nih-Image from Wayne Rasband); the following 640x480 bytes are the
c contents of the image.
c The total raw image size is 640*480 = 307200 bytes.
c The length of the Tiff header is 768 bytes, which is an entire fraction of 307200 bytes
c (768*400=307200).
c Therefore, in order to read the image by direct access the record length recl can be 768 bytes.
c However, because this is a very long word, I choose half of this length, that is recl=384.
c Then the number n of records needed to read the entire image will be twice as large, e.g. n=800.
c
integer n
character*20 list ! File list containing names of images to be renamed.
character*5 input ! Name of image file to be renamed xxxxr
character*9 output ! Renamed output: xxxx.raw
character*384 byte ! Variable for the record with length 384.
c
open(15,file=list,status='old')
5 read(15,6,end=99) input
6 format(a5) ! Five characters per input file name
c
c determine name of output file:
c
output=input//'.RAW'
c
c now open existing image files and files under new name output:
c
open(16,file=input,access='direct',
1 recl=384,status='old')
open(17,file=output,access='direct',
1 recl=384,status='new')
write(9,*) '. . .renaming raw files. . .'
c
c read contents of existing images and write them to output
c
do 10 n=1,800 ! n= number of record with recl=384.
read(16,'(a384)',rec=n) byte
write(17,'(a384)',rec=n) byte
10 continue
close(17)
close(16)
goto 5
c
99 write(9,*) 'All images renamed'
pause 100
100 continue
stop
end

6.15. Listing for program Rename_Tiff.f

program rename_Tiff
c
c Program to rename image files.
c By Michael Knappertsbusch, 11.8.2004.
c
c The program reads the contents of the files and writes them under a new name to disk.
c
c
c Explanations:
c Images have size of 640x480 pixels, 8 bit depth (256 grey levels), Tiff format.
c By standard the first 768 bytes of a Tiff file belong to the Tiff header (according
c to Documentation of Nih-Image from Wayne Rasband); the following 640x480 bytes are the
c contents of the image.
c The total image size is therefore 768+640*480 = 307968 bytes.
c The length of the Tiff header (768 bytes) is an entire fraction of 307968 bytes (768*401=307968).
c Therefore, in order to read the image by direct access the record length recl can be 768 bytes.
c However, because this is a very long word, I choose half of this length, that is recl=384.
c Then the number n of records needed to read the entire image will be twice as large, e.g. n=802.
c
integer n
character*20 list ! File list containing names of images to be renamed.
character*4 input ! Name of image file to be renamed
character*9 output ! Renamed input
character*384 byte ! Variable for the record with length 384.
c
open(15,file=list,status='old')
5 read(15,6,end=99) input
6 format(a4)
c
c determine name of output file:
c
output=input//'.TIFF'
c
c now open existing image files and files under new name output:
c
open(16,file=input,access='direct',
1 recl=384,status='old')
open(17,file=output,access='direct',
1 recl=384,status='new')
c
c read contents of existing images and write them to output
c
write(9,*) '. . .renaming tiff files. . .'
do 10 n=1,802 ! n= number of record with recl=384.
read(16,'(a384)',rec=n) byte
write(17,'(a384)',rec=n) byte
10 continue
close(17)
close(16)
goto 5
c
99 write(9,*) 'All images renamed'
pause 100
100 continue
stop
end

6.16. Listing for program Grid2.1.f

PROGRAM GRID
C
C Version 2.1, 6.4.1998, by M. Knappertsbusch.
C Modified 18. August 1999 from Version 2.0.
C
C . . .Batch operation mode. . .
C
C Verbesserte Version: Dimension der Matrizen erhoeht
C (kann nun bis 60x60 gridding files machen).
C DX und DY koennen jetzt variabel gewaehlt werden.
C Bereich des Gridding geht von 0 - XMax in X Richtung,
C und von 0 - YMax in Y Richtung.
C Ausgabe wurde verbessert.
C
C Modification 18. August 1999: Removed header info from
C "rotated files", program usage is now less confusing to
C other persons.
C
C Gridding program for C. leptoporus.
C Program calculates bivariate frequencies from X,Y
C data. Binwidth in X is DX (diameter), and in Y
C it is DY (elements).
C
C It applies the same grid as in Contour plot by
C Dave Lazarus.
C
C Arrays: FREQ(I,J), with I=Diameter (=X) and J=Elements (=Y).
C
C LABELX andLABELY are to indicate the categories in the
C output file.
C
C
INTEGER FREQ(1:60,1:60)
REAL LABELX(1:60),LABELY(1:60)
REAL XMAX,YMAX
INTEGER I,J,K
INTEGER NX,NY
REAL X,Y,DX,DY
CHARACTER*20 LIST
CHARACTER*24 INPUT
CHARACTER*25 STRING,OUTPUT
C
C REAL C1,C2,PHI
C CHARACTER*24 HEADER1
C CHARACTER*54 HEADER2
C CHARACTER*1 HEADER3
C
WRITE(9,*) '. . .Gridding X,Y data, version 2.0. . .'
WRITE(9,*) ' Batch operation mode'
WRITE(9,*) ' '
WRITE(9,*) 'The program reads a file containing the'
WRITe(9,*) 'names of all files to be gridded'
WRITE(9,*) '(The individual filenames must have names'
WRITE(9,*) 'WITH EXACTLY 20 CHARACTERS)'
WRITE(9,*) ' '
C
WRITE(9,*) 'Enter upper limit of X and Y (Xmax,Ymax)'
READ(9,*) XMAX,YMAX
WRITE(9,*) 'Enter intervals (DX,DY)'
READ(9,*) DX,DY
C
C Determine number of intervals in X and Y direction
C
NX=NINT(XMAX/DX)
NY=NINT(YMAX/DY)
WRITE(9,*) NX,' X-intervals'
WRITE(9,*) NY,' Y-intervals'
C
C File name handling:
C
WRITE(9,*) '. . .Enter list with sample names. . .'
WRITE(9,*) ' (Max. 20 chars)'
READ(9,1) LIST
1 FORMAT(A20)
C
C Open that file list and read sample names as variable STRING, and derive names of
C INPUT and OUTPUT files from STRING (K indicates the number of files gridded):
C
K=0
OPEN(14,FILE=LIST,STATUS='OLD')
2 READ(14,3,END=99) STRING
3 FORMAT(A25)
C
INPUT=STRING(1:20)
OUTPUT=STRING(1:20)//'.gridded'
C
C Now, open current input and output files:
C
OPEN(15,FILE=INPUT,STATUS='OLD')
OPEN(16,FILE=OUTPUT,STATUS='NEW')
C
C Reading header lines of rotated X,Y data:
C
C READ(15,6) HEADER1
C READ(15,7) HEADER2
C READ(15,*) C1,C2,PHI
C READ(15,8) HEADER3
C 6 FORMAT(A24)
C 7 FORMAT(A54)
8 FORMAT(A1)
C
C Initialize the gridding procedure:
C Initialize FREQ (I=Diameter, J=Elements):
C
DO 10,J=1,NY
DO 5, I=1,NX
FREQ(I,J)=0
5 CONTINUE
10 CONTINUE
C
C
C Determine Labels for X and Y values:
C
DO 11, I=1,NX
LABELX(I)=I*DX
11 CONTINUE
C
DO 12, I=1,NY
LABELY(I)=I*DY
12 CONTINUE
C
C
15 READ(15,*,END=999) X,Y
C
C Now gridding the data:
C I runs for Diameter, J runs for Elements.
C
DO 40, J=1,NY
DO 30, I=1,NX
IF (((J-1)*DY.LE.Y).AND.(Y.LT.J*DY)) THEN
IF (((I-1)*DX.LE.X).AND.(X.LT.I*DX)) THEN
FREQ(I,J)=FREQ(I,J)+1
END IF
END IF
30 CONTINUE
40 CONTINUE
GOTO 15
C
C Writing results to output file:
C
999 WRITE(16,50) OUTPUT
C WRITE(16,85) 'Center of rotation (X,Y),
C * rotation angle in degrees:'
C WRITE(16,*) C1,C2,PHI
WRITE(16,86) 'DeltaX, DeltaY, Number of X-intervals,
* Number Y-intervals:'
WRITE(16,*) DX,',',DY,',',NX,',',NY
WRITE(16,8) ' '
WRITE(16,70) ' ',(LABELX(I)-DX,'-',LABELX(I),I=1,NX)
WRITE(16,*) ' '
DO 60, J=1,NY
WRITE(16,80) LABELY(J)-DY,'-',LABELY(J),
*' : ',(FREQ(I,J),I=1,NX)
60 CONTINUE
C
CLOSE(15)
CLOSE(16)
K=K+1
GOTO 2
C
50 FORMAT(A25)
70 FORMAT(A3,60(F4.1,A1,F4.1,1X))
80 FORMAT(F4.1,A1,F4.1,A3,60(I4,1X))
85 FORMAT(A52)
86 FORMAT(A58)
C
99 WRITE(9,*) K,' files gridded'
STOP
END