' ************************************************************
' ** Name: CLAYES2K.BAS **
' ** Type: MAIN [ X ] FUNCTION [ ] SUB [ ] **
' ** Language: Microsoft QuickBASIC 4.00 or above **
' ** Library : Standard **
' ** Arguments: None - Main Module **
' ** Programmer: A. Eliason- Eliason Data- (508)477-1400 **
' ** Date: Version # in main module includes date **
' ** as: MMDDYYVV **
' ** Description: Please see comments within this module **
' ** **
' ** Usage: If no info here please see the line **
' ** comments. **
' ** **
' ** Y2K Check Year 2000 check performed 7-27-98 AHE **
' ** I have not found any Y2K Problems in this program **
' ************************************************************
' Clay Extrapolation Program - CLAYES2K.BAS
' Last rev 97-06-24 AHE Y2K 98-07-27
DEFINT I-N 'use Fortran conventions
DEFDBL A-H, O-Z
DIM PHIP(32), FP(32), FPcor(32), FPcor100(32), Iden(4), F(2002)
DIM XFit(4), YFit(4)
DIM s(32), index(32)
DIM Itext$(1000)
DIM Lable$(70)
DECLARE SUB CPI (PHIP(), FP(), N, DP, F(), IP, DPHI)
DECLARE SUB BPLOT (TMIN, TMAX, X(), N) 'EDS sub for FOPLOT
DECLARE SUB PressAnyKey ()
DECLARE SUB CSCOEF (N, X(), Y(), s(), index())
DECLARE SUB CSFIT1 (N, X(), Y(), s(), index(), X1, sx)
DECLARE SUB CSFIT2 (N, X1(), Y1(), YP1, YPN, X, Y)
DECLARE SUB FFIT (N, X(), Y(), IType, A, B)
DECLARE FUNCTION U2Phi (US)
Version$ = "072798.01 Y2K"
Diag% = 0
FileWrite% = 1
Cont% = 0
InterpFlag = 1
ISWFlag = 0 'no Schlee-Webster
Lprnt = 0
PhiNudge = .045 ' this factor is added to conversions from micron dia to Phi
' in order to adjust for the offset between whole micron
' diamaters and phi boundries eg: 1 micron actually= 9.965..Phi
' so we add .045 to Phi in order to cross the 10 Phi boundry.
NextFile:
CLS
PRINT
PRINT " USGS Sediment Size Analysis"
PRINT " Clay Fraction Estimation Version "; Version$
PRINT
PRINT "OPTIONS:"
INPUT " Do you want printed output (Y/N)? [N] ", A$
IF UCASE$(LEFT$(A$, 1)) = "Y" THEN Lprnt = 1
PRINT
ExType = 4 ' default to MEAN
InterpFlag = 1
' lets not kill this yet rev 6-18-97
'INPUT " Do you want to use Extrapolation (Y/N)? [Y] ", A$
'IF UCASE$(LEFT$(A$, 1)) = "N" THEN InterpFlag = 0
' If NO is selected then only the original data are printed/output
IF InterpFlag THEN
PRINT "Please READ the following carefully:"
PRINT
PRINT " You may elect to estimate clay content from the smallest measured size to"
PRINT " 13 phi. This may be done by LINEAR or EXPONENTIAL extrapolation, or the "
PRINT " MEAN of both sets of results. "
PRINT " Linear extrapolation may tend to over-estimate the clay present, while"
PRINT " exponential may under-estimate."
PRINT " It is up to the operator to determine the best method based"
PRINT " upon the original data."
PRINT
INPUT " Do you wish to use LINEAR, EXPONENTIAL, or MEAN (L/E/M)? [M] ", A$
IF UCASE$(LEFT$(A$, 1)) = "L" THEN ExType = 1
IF UCASE$(LEFT$(A$, 1)) = "E" THEN ExType = 3
' else ExType = the default, MEAN, = 4
PRINT
PRINT " The program needs to know the smallest measured size based upon"
PRINT " the Coulter callibration. Normally this will be a value between"
PRINT " 0.5 and 0.8 microns. You MUST enter a value."
Size:
PRINT
DO
INPUT " Smallest measured micron diameter"; UD
LOOP WHILE UD <= 0
IF UD < .5 OR UD > .8 THEN
PRINT USING " You have entered a value which corresponds to ###.## Phi"; U2Phi(UD)
INPUT " Is this correct (Y/N)? ", A$
IF UCASE$(LEFT$(A$, 1)) <> "Y" THEN GOTO Size:
END IF
PhiMin = U2Phi(UD) ' actulal phi size
PhiMin = PhiMin + PhiNudge ' correct to USGS boundry
IntPhi = INT(PhiMin) ' Integer portion of smallest measured
IntPhiFrac = IntPhi + 1 ' Fraction containing the smallest meas.
NdxPhiCor = IntPhiFrac + 6 ' index to Fraction Array
PhiCorFac = PhiMin - IntPhi ' divisor for correction
PctPhiCor = 100 - INT(PhiCorFac * 10000) / 100 ' just for information
PRINT USING " The value you entered corresponds to ##.## Phi"; PhiMin
PRINT USING " The ## Phi Fraction will be increased in order to account"; IntPhiFrac
PRINT USING " for the ##.## percent perceived to be undetected in this fraction."; PctPhiCor
IF Daig% THEN PRINT " Value will be divided by "; PhiCorFac
PRINT
END IF
PressAnyKey
' Read the input field lables
FOR I = 1 TO 68
READ Lable$(I)
IF Diag% = 1 THEN PRINT Lable$(I); " ~ ";
NEXT I
PRINT
' new file read routine 9-16-96 ############ AHE
FilePath:
INPUT "Enter Drive Letter for data file - [C:] ", Drv$
IF Drv$ = "" THEN Drv$ = "C:"
IF LEN(Drv$) = 1 THEN Drv$ = Drv$ + ":"
IF LEN(Drv$) > 2 THEN GOTO FilePath
PRINT
PRINT "You will now be asked for a PATH then a FILENAME. If you don't know, leave"
PRINT "one or both blank and you will be shown a Directory screen."
INPUT "File PATH to read? - eg '\DATA\' [\] ", FPN$
IF FPN$ = "" THEN FPN$ = "\"
IF LEFT$(FPN$, 1) <> "\" THEN FPN$ = "\" + FPN$
IF RIGHT$(FPN$, 1) <> "\" THEN FPN$ = FPN$ + "\"
PRINT "Data Path = "; Drv$ + FPN$
INPUT "FILENAME to read? "; FLN$
IF FLN$ = "" THEN ' Guess s/he doesn't know
FILES Drv$ + FPN$
PRINT "Please start again. When asked, enter a FILENAME from the list above."
PRINT "You MUST include the extension (eg AAAAAA.PG)" '6-97 AHE
PRINT "This program WILL NOT accept FILENAMEs of more than 6 characters nor ext-"
PRINT "tions greater than 3 characters"
GOTO FilePath
END IF
LabNo:
INPUT "Initial LAB NUMBER eg 'AHE001' ? ", LabNo$
LabNo$ = UCASE$(LabNo$)
LnLen = LEN(LabNo$)
FoundFlag% = 0
Drv$ = UCASE$(Drv$)
FPN$ = UCASE$(FPN$)
FLN$ = UCASE$(FLN$)
PRINT "Attempting to open "; Drv$ + FPN$ + FLN$
OPEN Drv$ + FPN$ + FLN$ FOR INPUT AS #1
DO WHILE NOT EOF(1)
INPUT #1, A$
IF (UCASE$(LEFT$(A$, LnLen)) = LabNo$) OR (UCASE$(MID$(A$, 2, LnLen)) = LabNo$) THEN
IF LEFT$(A$, 1) = "$" THEN A$ = MID$(A$, 2, LnLen)
FoundFlag% = 1
EXIT DO
END IF
IF Diag% THEN PRINT A$; " ";
LOOP
IF FoundFlag% THEN
PRINT
PRINT " Found "; A$
LN$ = A$
ELSE
BEEP
PRINT
PRINT " Did NOT find "; LabNo$
PRINT " Please try again with another Lab Number."
END IF
IF FoundFlag% = 0 THEN
CLOSE #1
PRINT
GOTO LabNo
END IF
INPUT "Do you wish to write extrapolated data to Disk? [Y]"; A$
IF UCASE$(LEFT$(A$, 1)) = "N" THEN FileWrite% = 0
IF FileWrite% = 1 THEN
RevType$ = "M_" 'default- mean of both
IF ExType = 1 THEN RevType$ = "L_" 'linear
IF ExType = 3 THEN RevType$ = "E_" 'exponential
OutFile$ = Drv$ + FPN$ + RevType$ + FLN$
IF LEN(FLN$) > 10 THEN
PRINT
PRINT "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!"
PRINT "CAN NOT Write output data to "; OutFile$
PRINT "!!!!!!!!!! File Name too Long !!!!!!!!!!"
PRINT "You must specify a new output file Path/Name!"
BEEP
ELSE
PRINT "Output data will be written to the following Path/File,"
PRINT "unless you enter another complete pathname."
END IF
PRINT "Output Filename [ "; OutFile$; " ] ";
INPUT "", A$
IF UCASE$(A$) = "Y" THEN A$ = ""
IF A$ <> "" THEN OutFile$ = A$
PRINT "Attempting to open "; OutFile$; " for output."
OPEN OutFile$ FOR OUTPUT AS #2
PRINT "Open sucessful."
' Write the attribute names to disk
FOR I = 1 TO 26
PRINT #2, Lable$(I); ",";
NEXT I
' attribute name # 27 is type of extrapolation
PRINT #2, "Extrapolation Method,";
' kill this LJP 10-24-96 PRINT #2, "Mode 1 Fq.%, Mode 2 Fq.%, Mode 3 Fq.%,Number of Modes,";
FOR I = 13 TO -5 STEP -1 ' -4 to -5 10/96 AHE
PRINT #2, I; " Phi Fq. %,";
NEXT I
PRINT #2, ""
END IF
PRINT
Cont% = 0 ' Correct error AHE 6-24-97
INPUT "Do you wish to pause and examine the data for each record? [Y]", A$
IF UCASE$(LEFT$(A$, 1)) = "N" THEN Cont% = 1 'True
'PressAnyKey
NextSamp:
Itext$(1) = LN$
I = 2
' Y2K WARNING - BELOW. Because input is alphanumeric, field # 8, year, should cause
' no problem, however it is formatted (2 digit 4 digit).
DO
A$ = "": B$ = ""
DO
DO ' ignore CR & LF
B$ = INPUT$(1, #1)
LOOP WHILE B$ = CHR$(10) OR B$ = CHR$(13)
IF B$ = "," THEN EXIT DO ' end of field - exit
A$ = A$ + B$ ' concatenate
IF B$ = "$" THEN EXIT DO ' because $ is not comma delimited
LOOP
Itext$(I) = A$
I = I + 1
IF Diag% THEN PRINT A$; "~";
LOOP UNTIL A$ = "$"
IF Diag% THEN
PRINT : PRINT
FOR I = 1 TO 68
PRINT USING "\ \&"; Lable$(I); Itext$(I)
IF (I MOD (22) = 0) THEN PressAnyKey
NEXT I
PressAnyKey
END IF
' check for proper field placement
IF Itext$(68) <> "$" THEN
PRINT : PRINT
PRINT "End of Record Character ($) NOT IN PROPER LOCATION"
PRINT "Aborting Program"
BEEP: BEEP: BEEP
PRINT : PRINT
END
END IF
' Now convert Cumulative Fq% to Fq%
FOR j = 1 TO 32
PHIP(j) = j - 6: FP(j) = 0: FPcor(j) = 0: FPcor100(j) = 0 'init the arrays
NEXT j
Cumul = 0: FqPct = 0: j = 1
FOR I = 67 TO 35 STEP -2 ' step from largest to smallest size
FqPct = VAL(Itext$(I)) - Cumul
FqPct = INT(FqPct * 100 + .499999999#) / 100' round to two places
Cumul = Cumul + FqPct
PHIP(j) = VAL(Itext$(I - 1)) ' Phi Size from input data
FP(j) = FqPct ' Actual phi percentage
IF Diag% THEN PRINT USING "## \ \ ##.## ###.##"; PHIP(j); Itext$(I); FP(j); Cumul
j = j + 1
'*** PHIP(j) now contains -5phi at 1, -4 at 2,... 11phi at 17 ***
NEXT I
PRINT
N = 20 ' -5 to 14 phi FRACTION inclusive
DPHI = .2
XSand = VAL(Itext$(17))
Grav = VAL(Itext$(18))
Silt = VAL(Itext$(19))
clay = VAL(Itext$(20))
' test for total percent = 100
CFP = 0
FOR j = 1 TO N
CFP = CFP + FP(j)
NEXT j
IF CFP > 101 OR CFP < 99 THEN
PRINT "Cumulative Frequency Percent is not 100 +/- 1"
PRINT
IF Lprnt THEN
LPRINT "Cumulative Frequency Percent is not 100 +/- 1"
LPRINT
END IF
PressAnyKey
GOTO NextSamp
END IF
' test if Phi intervals are even
DP = PHIP(2) - PHIP(1) '=.5
FOR j = 3 TO N
Test2 = PHIP(j) - PHIP(j - 1) '=.5
Test3 = ABS(DP - Test2) '=0
IF Diag% THEN PRINT DP, Test2, Test3: PRINT
IF (Test3 > .001) THEN
PRINT Iden$; " Phi intervals are not even"
PRINT
IF Lprnt THEN
PRINT Iden$; " Phi intervals are not even"
PRINT
END IF
PressAnyKey
GOTO NextSamp
END IF
DP = Test2
NEXT j
'PressAnyKey
RealStuff:
IF InterpFlag = 0 THEN GOTO PrintOut
PRINT
PRINT "Doing the initial Extrapolation"
FOR I = 1 TO N ' copy the array
FPcor(I) = FP(I)
NEXT I
FOR I = N TO 1 STEP -1 ' find the smallest non-zero fraction
IF FP(I) > 0 THEN
NZ = I ' NZ= index to smallest FRACTION with data
EXIT FOR
END IF
NEXT I
' Make the percentage correction to the fraction containg the smallest
' measurable data
imin = NdxPhiCor
IF FP(imin) > 0 THEN
FPcor(imin) = FP(imin) / PhiCorFac
PRINT USING "The smallest measurable fraction is ## Phi with a value of ##.##%"; IntPhiFrac; FP(imin)
PRINT USING "##.## has been changed to ##.##"; FP(imin); FPcor(imin)
PRINT
ELSE
PRINT USING "The ## Phi fraction was zero so it needed no change."; IntPhiFrac
PRINT
END IF
' Now do Extrapolation from the smallest non zero fraction to 14 FRACTION
' which is zero. Don't forget the 13 FRACTION contains the 12.001 to 13 data
PRINT USING "The smallest non-zero fraction is ## Phi with a value of ##.##%"; PHIP(NZ); FPcor(NZ)
IF ExType = 1 THEN PRINT "Linear ";
IF ExType = 3 THEN PRINT "Exponential ";
IF ExType = 4 THEN PRINT "Mean of L./E. ";
IF ExType <> 4 AND ExType <> 3 AND ExType <> 1 THEN PRINT "error at ExType. Program halted.": STOP
PRINT "Extrapolation will now be performed to 13 Phi.(14 Phi Fraction=0)"
PRINT
' Set up for calls to FFIT (Linear and Exponential)
index = NZ
XFit(1) = index
YFit(1) = FPcor(index)
XFit(2) = index
YFit(2) = FPcor(index)
XFit(3) = 20
YFit(3) = .0000001
CALL FFIT(3, XFit(), YFit(), 1, Alin, Blin)
CALL FFIT(3, XFit(), YFit(), 3, Aexp, Bexp)
index = NZ
'Npoints = 20 - NZ
'CLinFac = FPcor(index) / Npoints
FOR I = index + 1 TO 19
IF ExType = 1 THEN FPcor(I) = Alin + Blin * I ' linear
IF ExType = 3 THEN FPcor(I) = Aexp * EXP(Bexp * I) ' logarithmic
IF ExType = 4 THEN FPcor(I) = ((Alin + Blin * I) + (Aexp * EXP(Bexp * I))) / 2
NEXT I
' Re-normalize to 100 percent
TotalPct = 0
FOR I = 1 TO 20 ' N?
TotalPct = TotalPct + FPcor(I)
NEXT I
FOR I = 1 TO 20
FPcor100(I) = 100 * FPcor(I) / TotalPct
NEXT I
PrintOut:
PRINT
IF Cont% = 0 THEN PressAnyKey
'
' Print the header data
CLS
PRINT " USGS Sediment Size Analysis"
PRINT " Clay Fraction Estimation Version "; Version$
PRINT
PRINT "Lab # Field ID Proj. ID Cruise ID Requestor Latitude Longitude"
PRINT
PRINT USING "\ \"; Itext$(1); Itext$(2); Itext$(3); Itext$(4); Itext$(5); Itext$(9); Itext$(10)
PRINT
PRINT
PRINT "Device Location Depth Top Bottom Samp. Weight"
PRINT
PRINT USING "\ \"; Itext$(11); Itext$(12); Itext$(13); Itext$(14); Itext$(15); Itext$(16)
PRINT
IF Cont% = 0 THEN PressAnyKey
PRINT
PRINT " Original Revised";
IF ExType = 1 THEN PRINT " Linearly"
IF ExType = 3 THEN PRINT " Exponentially"
IF ExType = 4 THEN PRINT " by Mean of Linear and Exponential"
PRINT " Phi Size Freq. % Freq. %"
FOR j = 1 TO N
PRINT USING " ## ##.## ##.##"; PHIP(j); FP(j); FPcor100(j)
NEXT j
PRINT
IF Lprnt THEN
' Print the original data to LPT1
LPRINT " USGS Sediment Size Analysis"
LPRINT " Clay Fraction Estimation Version "; Version$; ""
LPRINT
LPRINT " Lab # Field ID Proj. ID Cruise ID Requestor Latitude Longitude"
LPRINT " ";
LPRINT USING "\ \"; Itext$(1); Itext$(2); Itext$(3); Itext$(4); Itext$(5); Itext$(9); Itext$(10)
LPRINT
LPRINT " Device Location Depth Top Bottom Samp. Weight"
LPRINT " ";
LPRINT USING "\ \"; Itext$(11); Itext$(12); Itext$(13); Itext$(14); Itext$(15); Itext$(16)
LPRINT
LPRINT " Original Revised";
IF ExType = 1 THEN LPRINT " Linearly"
IF ExType = 3 THEN LPRINT " Exponentially"
IF ExType = 4 THEN LPRINT " by Mean of Linear and Exponential"
LPRINT " Phi Size Freq. % Freq. %"
LPRINT
FOR j = 1 TO N
LPRINT USING " ## ##.## ##.##"; PHIP(j); FP(j); FPcor100(j)
NEXT j
LPRINT
END IF
' test for total percent = 100
CFP = 0
FOR j = 1 TO N
CFP = CFP + FP(j)
NEXT j
IF CFP > 101 OR CFP < 99 THEN
PRINT "Cumulative Frequency Percent is not 100 +/- 1"
PRINT
IF Lprnt THEN
LPRINT "Cumulative Frequency Percent is not 100 +/- 1"
LPRINT
END IF
PressAnyKey
GOTO NextSamp
END IF
' test if Phi intervals are even
DP = PHIP(2) - PHIP(1) '=.5
FOR j = 3 TO N
Test2 = PHIP(j) - PHIP(j - 1) '=.5
Test3 = ABS(DP - Test2) '=0
IF Diag% THEN PRINT DP, Test2, Test3: PRINT
IF (Test3 > .001) THEN
PRINT Iden$; " Phi intervals are not even"
PRINT
IF Lprnt THEN
PRINT Iden$; " Phi intervals are not even"
PRINT
END IF
PressAnyKey
GOTO NextSamp
END IF
DP = Test2
NEXT j
IF Cont% = 0 THEN PressAnyKey
Emp1 = PHIP(1) - .5 * DP 'Find center of 1st class interval
' to be used at end for statistics.
IP = (PHIP(N) - PHIP(1)) / DPHI + 1.5 'Calculate # of Extrapolated points
'(Span/Interp_Interval)+1.5 = # points
IF ISWFlag = 0 THEN
'PRINT "No S-W Extrapolation Requested"
'PRINT
'IF Lprnt THEN
' LPRINT "No S-W Extrapolation Requested"
' LPRINT
'END IF
IP = N ' Extrapolated Points = Original
DPHI = DP ' Delta Phi (interp) = Original
FOR I = 1 TO N
F(I) = FPcor100(I) ' Extrapolated Data = Original
'changed FP(i) to FPcor(i) 10-24-96 AHE
FP(I) = FPcor100(I) 'reload FP() 10-25-96
NEXT I
ELSE
PRINT "S-W Extrapolated Values will be used"
PRINT
IF Lprnt THEN
LPRINT "S-W Extrapolated Values were used"
LPRINT
END IF
FOR I = 1 TO N
FP(I) = FPcor100(I)
NEXT I
GOTO SchleeWebster
' Test CSCOEF ***************************************************
PRINT "Running CSCOEF ...";
CALL CSCOEF(N, PHIP(), FPcor(), s(), index())
PRINT " DONE"
FOR I = 1 TO N
PRINT s(I)
NEXT I
X1 = PHIP(1)
FOR I = 1 TO IP
'CALL CSFIT1(N, PHIP(), FPcor(), s(), index(), X1, SX)
'CALL CSFIT2(N, PHIP(), FPcor(), 0, 0, X1, SX)
F(I) = sx
X1 = X1 + DPHI
NEXT I
' Re-normalize to Phi Percentages
DPI = DPHI / DP
FOR I = 1 TO IP
IF F(I) < 0 THEN F(I) = 0
F(I) = F(I) * DPI
IF Diag% THEN
PRINT I, F(I)
IF (I MOD (15) = 0) THEN PressAnyKey
END IF
NEXT I
'Scale the data for plotting
Ymax = 0
FOR I = 1 TO IP
IF F(I) > Ymax THEN Ymax = F(I)
NEXT I
IF Lprnt THEN LPRINT "Plot Using CSFIT1 - Natural Cubic Spline"
'IF Lprnt THEN LPRINT "Plot Using CSFIT2 - Clamped Cubic Spline"
CALL BPLOT(0!, Ymax + 1.5, F(), IP)
IF Lprnt THEN LPRINT CHR$(12)
' Ends test ***********************************************************
SchleeWebster:
' move all the original values up by 1 in the arrays but leave first
' element as it was.
FOR I = 1 TO N
'IF Diag% THEN PRINT FP(J + 1), PHIP(J + 1)
j = N - I
FP(j + 2) = FP(j + 1)
PHIP(j + 2) = PHIP(j + 1)
NEXT I
PRINT
PRINT USING "### Points at ##.## Phi will be fit to ### points at #.### Phi"; N; DP; IP; DPHI
IF Lprnt THEN
LPRINT USING "### Points at ##.## Phi will be fit to ### points at #.### Phi"; N; DP; IP; DPHI
LPRINT
END IF
CALL CPI(PHIP(), FP(), N, DP, F(), IP, DPHI)
' Re-normalize to Phi Percentages
DPI = DPHI / DP
FOR I = 1 TO IP
IF F(I) < 0 THEN F(I) = 0
F(I) = F(I) * DPI
IF Diag% THEN
PRINT I, F(I)
IF (I MOD (15) = 0) THEN PressAnyKey
END IF
NEXT I
PRINT
END IF
'Here is a cut at the Lloyd Breslaus Sediment Class Program
'IF Cont% = 0 THEN PressAnyKey 'REM'd out to elininate double PAC 6-97
PRINT " Calculating Descriptive Classification Percentages"
PRINT
' Get the new percentages
Gravel = 0
Sand = 0
Silt = 0
clay = 0
FOR I = 1 TO 5 'Gravel
Gravel = Gravel + FP(I)
NEXT I
FOR I = 6 TO 10 'Sand
Sand = Sand + FP(I)
NEXT I
FOR I = 11 TO 14 'Silt
Silt = Silt + FP(I)
NEXT I
FOR I = 15 TO 20 'Clay
clay = clay + FP(I)
NEXT I
IF Sand = 0 THEN Sand = .001
IF Gravel = 0 THEN Gravel = .001
IF Silt = 0 THEN Silt = .001
IF clay = 0 THEN clay = .001
DO
IF Gravel >= 50 THEN C% = 1: EXIT DO
IF Gravel >= 10 THEN C% = 2: EXIT DO
XSand = Sand
Sand = Sand + Gravel
IF Sand >= 75 THEN C% = 3: EXIT DO
IF Silt >= 75 THEN C% = 4: EXIT DO
IF clay >= 75 THEN C% = 5: EXIT DO
' Sand, Silt, & Clay are all less than 75%
SanSil = Sand / Silt
ClySnd = clay / Sand
SilCly = Silt / clay
IF Sand <= 20 THEN
IF SanSil > 1 THEN C% = 6: EXIT DO
IF SilCly < 1 THEN C% = 7: EXIT DO
IF ClySnd > 1 THEN
C% = 8: EXIT DO
ELSE C% = 9: EXIT DO
END IF
ELSE ' Sand> 20%
IF clay <= 20 THEN
IF SanSil < 1 THEN C% = 9: EXIT DO
IF SilCly > 1 THEN C% = 10: EXIT DO
IF ClySnd > 1 THEN
C% = 6: EXIT DO
ELSE
C% = 11: EXIT DO
END IF
ELSE ' Clay > 20%
IF Silt > 20 THEN C% = 12: EXIT DO
IF ClySnd > 1 THEN
C% = 6: EXIT DO
ELSE
C% = 11: EXIT DO
END IF
END IF
END IF
PRINT "ERROR IN BRESLAU ROUTINE": BEEP: END ' ya can't get heah from theya
LOOP
SELECT CASE C%
CASE 1
C$ = "GRAVEL"
CASE 2
C$ = "GRAVELLY SEDIMENT"
CASE 3
C$ = "SAND"
CASE 4
C$ = "SILT"
CASE 5
C$ = "CLAY"
CASE 6
C$ = "SANDY CLAY"
CASE 7
C$ = "SILTY CLAY"
CASE 8
C$ = "CLAYEY SILT"
CASE 9
C$ = "SANDY SILT"
CASE 10
C$ = "SILTY SAND"
CASE 11
C$ = "CLAYEY SAND"
CASE 12
C$ = "SAND SILT CLAY"
CASE ELSE
C$ = "Error in classification routine"
END SELECT
PRINT " Original Revised"
PRINT USING "Gravel % ##.## ##.##"; VAL(Itext$(18)); Gravel
PRINT USING "Sand % ##.## ##.##"; VAL(Itext$(17)); XSand
PRINT USING "Silt % ##.## ##.##"; VAL(Itext$(19)); Silt
PRINT USING "Clay % ##.## ##.##"; VAL(Itext$(20)); clay
Total = Gravel + XSand + Silt + clay
PRINT USING " Total = ###.##"; Total
PRINT
PRINT "Revised Classification: "; C$
Verbal$ = C$
PRINT
IF Lprnt THEN
LPRINT " Original Revised"
LPRINT USING " Gravel % ##.## ##.##"; VAL(Itext$(18)); Gravel
LPRINT USING " Sand % ##.## ##.##"; VAL(Itext$(17)); XSand
LPRINT USING " Silt % ##.## ##.##"; VAL(Itext$(19)); Silt
LPRINT USING " Clay % ##.## ##.##"; VAL(Itext$(20)); clay
LPRINT USING " Total = ###.##"; Total
LPRINT
LPRINT " Revised Classification: "; C$
LPRINT
END IF
IF Cont% = 0 THEN PressAnyKey
'Scale the data for plotting
Ymax = 0
FOR I = 1 TO IP
IF F(I) > Ymax THEN Ymax = F(I)
NEXT I
CALL BPLOT(0!, Ymax + 1.5, F(), IP)
STATS:
' now do the statistics via Schlee/Webster 1966
' Start with MODE
'IP = N ' need if not using schlee-webster
Del1 = -.1
PRINT "Revised Statistics"
IF Lprnt THEN
LPRINT
LPRINT " Revised Statistics"
END IF
FOR I = 2 TO IP
Del2 = F(I) - F(I - 1)
IF ((Del2 * Del1) < 0) AND (Del2 < 0) AND ((F(I - 1) - 5! * DPHI) > 0) THEN
Xi = I - 2
Xem = Emp1 + Xi * DPHI
Xf = F(I - 1)
I = I + 1
PRINT " Mode "; Xem, Xf
IF Lprnt THEN
'LPRINT
'LPRINT " Mode "; Xem, Xf
END IF
END IF
Del1 = Del2
NEXT I
' MEDIAN
Sum = 0
FOR I = 1 TO IP
Sum = F(I) + Sum
IF (Sum - 50!) >= 0 THEN EXIT FOR
NEXT I
Xi = I - 1
Ct = Emp1 + Xi * DPHI
CtMed = Ct - (Sum - 50!) * DPHI / F(I) + .5 * DPHI
Fmt$ = " \ \ =#####.##"
PRINT USING Fmt$; "Median"; CtMed
IF Lprnt THEN LPRINT USING Fmt$; "Median"; CtMed
' Moments about the Mean
Sum = 0: Sum1 = 0: Sum2 = 0: Sum3 = 0: Sum4 = 0
FOR I = 1 TO IP
Xi = I - 1
Ct = Emp1 + Xi * DPHI
Sum = Sum + F(I)
Sum1 = Sum1 + F(I) * Ct
Sum2 = Sum2 + F(I) * Ct ^ 2
Sum3 = Sum3 + F(I) * Ct ^ 3
Sum4 = Sum4 + F(I) * Ct ^ 4
NEXT I
En1 = Sum1 / Sum
En2 = Sum2 / Sum
En3 = Sum3 / Sum
En4 = Sum4 / Sum
Zm2 = En2 - En1 ^ 2
Zm3 = En3 - 3! * En2 * En1 + 2! * En1 ^ 3
Zm4 = En4 + En1 * (-4 * En3 + 6! * En1 * En2 - 3! * En1 ^ 3)
DPHI2 = DPHI * DPHI
Zm4 = Zm4 - .5 * DPHI2 * Zm2 + .02916667# * DPHI2 * DPHI2
Zm2 = Zm2 - DPHI2 / 12!
'PRINT En1; En2; En3; En4; Zm2; Zm3; Zm4
' Standard Deviation
Sigma = SQR(Zm2)
' Skewness
Skew = .5 * Zm3 / (Sigma * Zm2)
' Kurtosis
ZKurt = Zm4 / Zm2 ^ 2 - 3!
' Print those puppies
PRINT USING Fmt$; "Mean"; En1
PRINT USING Fmt$; "Standard Deviation"; Sigma
PRINT USING Fmt$; "Skewness"; Skew
PRINT USING Fmt$; "Kurtosis"; ZKurt
PRINT
IF Lprnt THEN
LPRINT USING Fmt$; "Mean"; En1
LPRINT USING Fmt$; "Standard Deviation"; Sigma
LPRINT USING Fmt$; "Skewness"; Skew
LPRINT USING Fmt$; "Kurtosis"; ZKurt
LPRINT CHR$(12) ' Form Feed
END IF
' write the corrected data to disk
IF FileWrite% THEN
FOR I = 1 TO 16
PRINT #2, Itext$(I); ","; ' Lab# - Samp Wt.
NEXT I
PRINT #2, USING "###.##,"; Sand; Gravel; Silt; clay;
PRINT #2, Verbal$; ",";
PRINT #2, USING "###.##,"; CtMed; En1; Sigma; Skew; ZKurt;
'Print #2, "9999,9999,9999,0," ' NEEDS WORK AHE - Killed LJP 10/96
IF ExType = 1 THEN PRINT #2, "LINEAR,";
IF ExType = 3 THEN PRINT #2, "EXPONENTIAL,";
IF ExType = 4 THEN PRINT #2, "MEAN_L_E,";
FOR I = 19 TO 1 STEP -1
PRINT #2, USING "###.##,"; FPcor100(I); 'round to 2 dp'
NEXT I
PRINT #2, "$"
END IF
IF Cont% = 0 THEN PRINT "For next Sample - "; : PressAnyKey
' See if the file is at end of data
DO WHILE NOT EOF(1)
INPUT #1, LN$
IF LN$ <> "" THEN GOTO NextSamp ' start next record
LOOP
CLS
CLOSE #1
CLOSE #2
PRINT
PRINT "End of File"
INPUT "Enter 'C' to Continue or 'Q' to quit program. [Q] ", A$
IF UCASE$(LEFT$(A$, 1)) = "C" THEN RESTORE: GOTO NextFile
END
' Label Data
' Y2K WARNING - Need to check if Year will crash if four digit or zero
DATA Lab Number,Field ID,Project ID,Cruise ID,Requestor,Month,Day,Year
DATA Latitude,Longitude,Sample Device,Location,Depth,Top,Bottom,Sample Weight
DATA % Sand,% Gravel,% Silt,% Clay,Verbal Equivalent,Median,Mean
DATA Standard Deviation,Skewness,Kurtosis,Mode 1,Mode 1 Frequency %
DATA Mode 2,Mode 2 Frequency %,Mode 3,Mode 3 Frequency %,Number of Modes
DATA "11 Phi","11 Phi Cumulative Fq %","10 Phi","10 Phi Cummulative Fq %"
DATA " 9 Phi"," 9 Phi Cumulative Fq %"," 8 Phi"," 8 Phi Cummulative Fq %"
DATA " 7 Phi"," 7 Phi Cumulative Fq %"," 6 Phi"," 6 Phi Cummulative Fq %"
DATA " 5 Phi"," 5 Phi Cumulative Fq %"," 4 Phi"," 4 Phi Cummulative Fq %"
DATA " 3 Phi"," 3 Phi Cumulative Fq %"," 2 Phi"," 2 Phi Cummulative Fq %"
DATA " 1 Phi"," 1 Phi Cumulative Fq %"," 0 Phi"," 0 Phi Cummulative Fq %"
DATA "-1 Phi","-1 Phi Cumulative Fq %","-2 Phi","-2 Phi Cummulative Fq %"
DATA "-3 Phi","-3 Phi Cumulative Fq %","-4 Phi","-4 Phi Cummulative Fq %"
DATA "-5 Phi","-5 Phi Cumulative Fq %",End Flag
DATA ""
' ************************************************************
' ** Name: BPLOT **
' ** Type: MAIN [ ] FUNCTION [ ] SUB [ X ] ** **
' ** Language: Microsoft QuickBASIC 4.00 or above **
' ** Library : Standard **
' ** Arguments: **
' ** Programmer: A. Eliason- Eliason Data- (508)477-1400 **
' ** Date: Version # in main module includes date **
' ** as: MMDDYYVV **
' ** Description: Please see comments within this module **
' ** **
' ** Usage: If no info here please see the line **
' ** comments. **
' ************************************************************
'
'
'
'
'
' ************************************************************
' ** Name: BPLOT **
' ** Type: MAIN [ ] FUNCTION [ ] SUB [ X ] ** **
' ** Language: Microsoft QuickBASIC 4.00 or above **
' ** Programmer: A. Eliason- Eliason Data- (508)477-1400 **
' ** Date: Version # in main module includes date **
' ** as: MMDDYYVV **
' ** Description: Please see comments within this routine **
' ** **
' ** Usage: If no info here please see the line **
' ** comments. **
' ************************************************************
SUB BPLOT (TMIN, TMAX, X(), N)
SHARED Lprnt, Cont%
M = 0
SPAN = TMAX - TMIN
UNIT = SPAN / 70
FOR I = 1 TO N
M = M + 1
IF M = 20 THEN M = 0: IF Cont% = 0 THEN PressAnyKey
ISPC = (X(I) - TMIN) / UNIT
Fract = (I - 1) / 5 - 5
IFract = Fract - .4999
'PRINT USING "## ##.##### "; IFract; X(i);
PRINT USING "##.##### "; X(I);
PRINT SPC(ISPC); "*" '; ' SPC(70 - ISPC);
'IF Lprnt THEN
' LPRINT USING "## ##.##### "; IFract; X(i);
' LPRINT SPC(ISPC); "*" '; ' SPC(70 - ISPC);
'END IF
NEXT I
PRINT
END SUB
SUB CPI (H(), Q(), N, DELH, F(), M, DELHC)
SHARED Diag%
Q(1) = 0
Q(N + 2) = 0
FOR j = 2 TO N
DELQ = (Q(j + 1) - Q(j)) / DELH
Y = 3! * Q(j) - 3! * Q(j + 1) + Q(j + 2)
Z = H(j)
IMAX = DELH / DELHC '= number of Extrapolated points per interval
FOR I = 1 TO IMAX
JP = (Z - H(2)) / DELHC + 1.4999 ' Prevent rounding error
'PRINT "*********** JP= "; JP; " ": INPUT A$
XK = (Z - H(j)) / DELH
X = XK + 1!
A = Q(j - 1) + XK * (Y - Q(j - 1))
FC = A + X * (Q(j) - A + .5 * (Q(j + 1) - 2! * Q(j) + A) * XK)
IF (j - 2) = 0 THEN GOTO Cont12
IF (j - N) <> 0 THEN GOTO Cont5
Cont12:
' Test for Extrapolated value less than zero
'PRINT FC
IF FC < 0 THEN
Z = H(j) + DELHC
IF (j - 2) = 0 THEN
Q(j) = .01
PRINT "Exponential Extrapolation at Beginning"
GOTO Cont7
END IF
IF (j - N) = 0 THEN
Q(j + 1) = .01
PRINT "Exponential Extrapolation at End"
ELSE
GOTO Cont5
END IF
Cont7:
IF Q(j + 1) = 0 THEN Q(j + 1) = .000001 'prevent div by zero
IF Q(j) = 0 THEN Q(j) = .000001
SLPE = LOG(Q(j) / Q(j + 1)) / DELH
FOR II = 2 TO IMAX
JP = (Z - H(2)) / DELHC + 1.4999
ANTI = SLPE * (Z - H(j))
F(JP) = Q(j) / EXP(ANTI)
Z = Z + DELHC
NEXT II
GOTO Cont1:
END IF
Cont5: F(JP) = FC
Cont2: Z = Z + DELHC
NEXT I
Cont1:
NEXT j
F(JP + 1) = 0
IF Diag% THEN PRINT "CPI WAS RUN. EXITING CPI."
END SUB
SUB CSCOEF (N, X(), Y(), s(), index()) STATIC
' Copyright (c) 1993 Crescent Software
' Cublic spline coefficients subroutine
' Input
' N = number of X and Y data points
' X() = array of X data points (N rows by 1 column)
' Y() = array of Y data points (N rows by 1 column)
' Output
' S() = array of cubic spline coefficients
' (N rows by 1 column)
' INDEX() = array of indices (N rows by 1 column)
REDIM RHO(N), TAU(N)
NM1 = N - 1
FOR I = 1 TO N
index(I) = I
NEXT I
' ascending order data sort
FOR I = 1 TO NM1
IP1 = I + 1
FOR j = IP1 TO N
II = index(I)
IJ = index(j)
IF (X(II) > X(IJ)) THEN
ITEMP = index(I)
index(I) = index(j)
index(j) = ITEMP
END IF
NEXT j
NEXT I
NM2 = N - 2
RHO(2) = 0#
TAU(2) = 0#
FOR I = 2 TO NM1
IIM1 = index(I - 1)
II = index(I)
IIP1 = index(I + 1)
HIM1 = X(II) - X(IIM1)
HI = X(IIP1) - X(II)
Temp = (HIM1 / HI) * (RHO(I) + 2#) + 2#
RHO(I + 1) = -1# / Temp
D = 6# * ((Y(IIP1) - Y(II)) / HI - (Y(II) - Y(IIM1)) / HIM1) / HI
TAU(I + 1) = (D - HIM1 * TAU(I) / HI) / Temp
NEXT I
s(1) = 0#
s(N) = 0#
' compute cubic spline coefficients
FOR I = 1 TO NM2
IB = N - I
s(IB) = RHO(IB + 1) * s(IB + 1) + TAU(IB + 1)
NEXT I
ERASE RHO, TAU
END SUB
SUB CSFIT1 (N, X(), Y(), s(), index(), X1, sx) STATIC
' Copyright (c) 1993 Crescent Software
' Cublic spline Extrapolation subroutine
' Input
' N = number of X and Y data points
' X() = array of X data points (N rows by 1 column)
' Y() = array of Y data points (N rows by 1 column)
' S() = array of cubic spline coefficients
' (N rows by 1 column)
' INDEX() = array of indices (N rows by 1 column)
' X1 = X data value to fit
' Output
' SX = cubic spline Extrapolated value for X
FOR I = 2 TO N
II = index(I)
IF (X1 <= X(II)) THEN EXIT FOR
NEXT I
L = I - 1
IL = index(L)
ILP1 = index(L + 1)
A = X(ILP1) - X1
B = X1 - X(IL)
HL = X(ILP1) - X(IL)
sx = A * s(L) * (A * A / HL - HL) / 6# + B * s(L + 1) * (B * B / HL - HL) / 6# + (A * Y(IL) + B * Y(ILP1)) / HL
END SUB
SUB CSFIT2 (N, X1(), Y1(), YP1, YPN, X, Y) STATIC
' Copyright (c) 1989, 1990 Crescent Software
' Clamped cubic spline Extrapolation subroutine
' Input
' N = number of X and Y data points
' X1() = vector of X data points ( N rows )
' Y1() = vector of Y data points ( N rows )
' YP1 = derivative at data point 1
' YPN = derivative at data point N
' X = X data point to fit
' Output
' Y = Extrapolated Y data point
' NOTE: X1(1) < X1(2) < X1(3) < ... < X1(N)
REDIM U(N), YPP(N)
YPP(1) = -.5
U(1) = (3# / (X1(2) - X1(1))) * ((Y1(2) - Y1(1)) / (X1(2) - X1(1)) - YP1)
FOR I = 2 TO N - 1
SIG = (X1(I) - X1(I - 1)) / (X1(I + 1) - X1(I - 1))
P = SIG * YPP(I - 1) + 2#
YPP(I) = (SIG - 1#) / P
U(I) = (6# * ((Y1(I + 1) - Y1(I)) / (X1(I + 1) - X1(I)) - (Y1(I) - Y1(I - 1)) / (X1(I) - X1(I - 1))) / (X1(I + 1) - X1(I - 1)) - SIG * U(I - 1)) / P
NEXT I
QN = .5
UN = (3# / (X1(N) - X1(N - 1))) * (YPN - (Y1(N) - Y1(N - 1)) / (X1(N) - X1(N - 1)))
YPP(N) = (UN - QN * U(N - 1)) / (QN * YPP(N - 1) + 1#)
FOR K = N - 1 TO 1 STEP -1
YPP(K) = YPP(K) * YPP(K + 1) + U(K)
NEXT K
' Extrapolate
KLO = 1
KHI = N
WHILE (KHI - KLO > 1)
K = (KHI + KLO) / 2#
IF (X1(K) > X) THEN
KHI = K
ELSE
KLO = K
END IF
WEND
H = X1(KHI) - X1(KLO)
A = (X1(KHI) - X) / H
B = (X - X1(KLO)) / H
Y = A * Y1(KLO) + B * Y1(KHI) + ((A ^ 3 - A) * YPP(KLO) + (B ^ 3 - B) * YPP(KHI)) * (H ^ 2) / 6#
ERASE U, YPP
END SUB
SUB FFIT (N, X(), Y(), IType, A, B) STATIC
' Function curve fit subroutine
' Input
' N = number of data points (N >= 3)
' X() = vector of X data points ( N rows )
' Y() = vector of Y data points ( N rows )
' ITYPE = type of curve fit
' 1 = linear Y = A + B * X
' 2 = logarithmic Y = A + B * LOG(X)
' 3 = exponential Y = A * EXP(B * X)
' Output
' A, B = coefficients of curve fit
X1 = 0#
Y1 = 0#
Z = 0#
X2 = 0#
Y2 = 0#
SELECT CASE IType
CASE 1
' Linear
FOR I = 1 TO N
X1 = X1 + X(I)
Y1 = Y1 + Y(I)
X2 = X2 + X(I) * X(I)
Y2 = Y2 + Y(I) * Y(I)
Z = Z + X(I) * Y(I)
NEXT I
X1 = X1 / N
Y1 = Y1 / N
B = (Z - N * X1 * Y1) / (X2 - N * X1 * X1)
A = Y1 - B * X1
CASE 2
' Logarithmic
FOR I = 1 TO N
XLOGX = LOG(X(I))
X1 = X1 + XLOGX
Y1 = Y1 + Y(I)
X2 = X2 + XLOGX * XLOGX
Y2 = Y2 + Y(I) * Y(I)
Z = Z + XLOGX * Y(I)
NEXT I
X1 = X1 / N
Y1 = Y1 / N
B = (Z - N * X1 * Y1) / (X2 - N * X1 * X1)
A = Y1 - B * X1
CASE 3
' Exponential
FOR I = 1 TO N
XLOGY = LOG(Y(I))
X1 = X1 + X(I)
X2 = X2 + X(I) * X(I)
Y1 = Y1 + XLOGY
Y2 = Y2 + XLOGY * XLOGY
Z = Z + X(I) * XLOGY
NEXT I
X1 = X1 / N
Y1 = Y1 / N
B = (Z - N * X1 * Y1) / (X2 - N * X1 * X1)
A = EXP(Y1 - B * X1)
END SELECT
END SUB
DEFSNG A-H, O-Z
SUB PressAnyKey
X% = POS(0)
Y% = CSRLIN
PRINT "Press any key to continue.";
DO
LOOP WHILE INKEY$ = ""
LOCATE Y%, X%
PRINT " ";
LOCATE Y%, X%
END SUB
DEFDBL A-H, O-Z
FUNCTION U2Phi (Usize)
' convert micron size to phi
phi = -2.12720466241656D-16 + (-1.44269504088896#) * LOG(Usize / 1000)
'PRINT USING "###.######### mm = "; Usize;
'PRINT USING "###.### phi"; Phi
U2Phi = phi
END FUNCTION