Subroutine Idget(Ix,Iy,Iopt,Ifit,Iorder,Npt,Xvals,Numbr,Data, $ Iflag,Xloc,Nxloc) c c *********************************************************************** c * * c * Subroutine IDGET - To get from the user the X and Y columns * c * to integrate or differentiate, the points to evaluate it * c * at, and the curve fit to evaluate it with. * c * * c * Output values: * c * * c * Ifit - The curve fit for integration/differentiation * c * Iflag - Error flag to indicate the user wishes to abort * c * Iopt - Flag to indicate whether we integrate or * c * differentiate * c * Iorder - The degree of the least squares fit if chosen * c * Ix - The column number of the X values to operate on * c * Iy - The column number of the Y values to operate on * c * Xvals - The X locations to evaluate the integral/derivative * c * at * c * X - The X values to operate on (returned in common * c * block MAIN * c * Y - The Y values to operate on (returned in common * c * block MAIN * c * * c * Programmed by Steven C. Fisher, DTNSRDC Code 1522 * c * * c *********************************************************************** c Include 'Cmain' c Virtual Data(1) c Dimension Xvals(1),Xloc(1) c c Initialize error flag (no error = 0) c Iflag = 0 c c See if there is any data in the program c If (Ncols .EQ. 0) Then Write(7,100) 100 Format(' ERROR - No Data in the program - routine aborted') Iflag = 1 Return End If c 10 Continue c c Input the curve to integrate/differentiate c Write(7,150) 27,27,27 150 Format(1X,1A1,'[H',1A1,'[J'/1X,1A1,'#6 DATA INTEGRATION/', $ 'DIFFERENTIATION ROUTINE'/) 15 Continue Write(7,200) 200 Format(' Enter the column numbers corresponding to', $ ' the X and Y values'/) Call Indata(2,'II',Ix,Iy) c c See if the column numbers are within the limits, and c if both columns have the same number of points. c Call Colchk(Ix,Iy,Ifl,7) If (Ifl .EQ. 1) Go To 15 c c See if only one point is in each column (error condition) c If (Npoint(Ix) .EQ. 1) Then Write(7,250) 250 Format(' Only one point in each column - routine aborted') Iflag = 1 Return End If c c See if the user likes their column choices c Write(7,300) 300 Format(' Are these column numbers correct? (y/n) '$) Call Instrg('YyNn',4,Ifl) If (Ifl .GT. 2) Then Write(7,350) 350 Format(' Do you wish to stop? (y/n) '$) Call Instrg('YyNn',4,Ifl) If (Ifl .GT. 2) Go To 10 Iflag = 1 Return End If c c Determine if the user will integrate or differentiate c Write(7,400) 400 Format(' Select the calculus operation:'//' 1) Integration'/ $ ' 2) First Derivative'// $ ' Enter 1-2 ',$) Call Inchek(1,2,Iopt) c c Determine the way the X values to integrate or c differentiate at will be input c Write(7,450) 450 Format(' X Location Selection....Enter Either'// $ ' 1) X values stored in the given X column'/ $ ' 2) Specific X values'/ $ ' 3) Minimum, maximum, and increment values'/ $ ' 4) X values stored in a different column'// $ ' Enter 1-4 '$) Call Inchek(1,4,Ifl) c c Enter the values to integrate/differentiate c If (Ifl .EQ. 1) Then Call Getvec(Data,Npoint,Ix,Xvals,Numbr) Else If (Ifl .EQ. 2) Then c c Specific X values c Write(7,500) 500 Format(' Enter the X values - enter "S" to stop') Call INPT(Xvals,Numbr,If,1) Else If (Ifl .EQ. 3) Then c c Minumum, maximum, and increment c 20 Continue Write(7,550) 550 Format(' Enter the minimum, maximum, and increment X values') Call Indata(3,'RRR',Xmin,Xmax,Xinc) c Read(5,*,END=20,ERR=20) Xmin,Xmax,Xinc c c Compute the X values c Numbr = 1 + (Xmax - Xmin) / Xinc If (Numbr .GT. Nptmax) Numbr = Nptmax Do 25 I=1,Numbr Xvals(I) = Xmin + (I-1) * Xinc 25 Continue Else c c X values stored in a column; get the column number from c the user; diaplay column title, and have user ok it. c 30 Continue Write(7,600) 600 Format(' Enter the column number containing the X values to', $ ' evaluate at ',$) Call Inchek(1,Ncols,Icol) Call Name(Titl,Icol,Title,Ttllen) Write(7,650) (Titl(I),I=1,Ttllen(Icol)) 650 Format(/' Column title is : ',A1/' Is it', $ ' correct? (y/n) ',$) Call Instrg('YyNn',4,Iflagg) If (Iflagg .GT. 2) Go To 30 Call Getvec(Data,Npoint,Icol,Xvals,Numbr) End If c Kflag = 0 Do 32 I=2,Numbr If (Xvals(I) .LT. Xvals(I-1)) Kflag = 1 32 Continue c If (Kflag .GT. 0 .AND. Iopt .EQ. 2) Then Write(7,690) 690 Format(' Sort the X values to interpolate at? (y/n) ',$) Call Instrg('YyNn',4,I) If (I .LT. 3) Call Sort(Xvals,Y,Numbr) Else If (Iopt .EQ. 1) Then Call Sort(Xvals,Y,Numbr) End If c c Get the data curve fit type c Write(7,700) 700 Format(' Enter the curve fit to integrate/differentiate:'// $ ' 1) Straight line'/ $ ' 2) Least squares polynomial'/ $ ' 3) Cubic spline'/ $ ' 4) Parametric spline - type A'/ $ ' 5) Parametric spline - type B'/ $ ' 6) Parametric spline - type C'/ $ ' 7) Least squares cubic spline (type A)'/ $ ' 8) Least squares cubic spline (type B)'// $ ' Enter 1 - 8 '$) Call Inchek(1,8,Ifit) c c Get the X and Y values stored in the columns c Call Getvec(Data,Npoint,Ix,X,Npt) Call Getvec(Data,Npoint,Iy,Y,Npt) c c Sort the data (unless a parametric spline of type C c is selected) c If (Ifit .NE. 6) Then Call Sort(X,Y,Npt) End If c c If a least squares fit is selected, get the degree of the polynomial c If (Ifit .EQ. 2) Then c Write(7,750) 750 Format(' Enter the degree of the polynomial (3 = cubic) '$) Call Inchek(1,20,Iorder) End If c c Least squares cubic spline fit c If (Ifit .GE. 7) Then 35 Continue Write(7,800) 800 Format(/' Least Sqaures Cubic Spline Fit:'/ $ ' Enter the X values of the junction points and',/ $ ' end points of the spline segments'// $ ' Select:'// $ ' 1) X values input by keyboard'/ $ ' 2) X values stored in a column'/ $ ' 3) Stop'// $ ' Enter 1-3: ',$) Call Inchek(1,3,Ich) c c Get the segment junction points c If (Ich .EQ. 1) Then Write(7,500) Call INPT(Xloc,Nxloc,If,1) Else If (Ich .EQ. 2) Then 40 Continue Write(7,850) 850 Format(' Enter the column number containing the X juncture', $ ' values ',$) Call Inchek(1,Ncols,Icol) Call Name(Titl,Icol,Title,Ttllen) Write(7,650) (Titl(I),I=1,Ttllen(Icol)) Call Instrg('YyNn',4,I) If (I .GT. 2) Go To 40 Call Getvec(Data,Npoint,Icol,Xloc,Nxloc) End If c If (Ich .EQ. 3 .OR. Nxloc .LT. 2) Then Iflag = 1 Return End If c If (Nxloc .GT. 60) Then Write(7,900) 900 Format(' ERROR - Too many points in X juncture column -', $ ' Try Again') Go To 35 End If c c Make certain that there are datapoints between the junction points c Call Sort(Xloc,Y,Nxloc) Call Getvec(Data,Npoint,Ix,X,Npt) Call Getvec(Data,Npoint,Iy,Y,Npt) Call Sort(X,Y,Npt) c Do 60 I=1,Nxloc-1 Istart = -1 Iend = Npt Do 50 J=1,Npt If (Xloc(I) .LE. X(J) .AND. Istart .EQ. -1) Istart = J If (Xloc(I+1) .LT. X(J)) Then Iend = J - 1 Go To 55 End If 50 Continue 55 Continue Istart = Max0(1,Istart) Ntot = Iend - Istart + 1 If (Ntot .EQ. 0) Then Write(7,950) Xloc(I),Xloc(I+1) 950 Format(' ERROR - There are no datapoints between the'/ $ ' junction points X= ',G12.5,' and X= ',G12.5/) c Go To 35 End If 60 Continue End If c Return End Subroutine Tabin(Iflag) c c *********************************************************************** c * * c * Subroutine Tabin - To input the column numbers and their * c * corresponding header labels. * c * * c * Programmed by Steven C. Fisher, DTNSRDC Code 1522 * c * * c *********************************************************************** c Include 'Cmain' Include 'Ctable' c Call Point(Loop,Ntab+1,Tcolen) Tcolen(Ntab+1) = 0 c c If this isn't the title for the table, write out c instructions to the user c c If (Ntab .GT. 0) Then c 10 Continue c c Get the column number for this table column c Write(7,100) Ntab 100 Format(/' Enter the column number for table column',I3,' (', $ 'enter "s" to stop) ',$) Call Inpt(X,N,Ifl,2) c c Check the input column number for errors c If (N .EQ. 0) Then Iflag = 1 Return End If c Tabcol(Ntab+1) = X(1) c If (Tabcol(Ntab+1) .LT. 1 .OR. Tabcol(Ntab+1) .GT. Ncols) Then Write(7,200) 200 Format(' ERROR - Column number is out of range - Try Again') Go To 10 End If c c See if the user likes their column choices c Call Name(Titl,Tabcol(Ntab+1),Title,Ttllen) Write(7,250) (Titl(I),I=1,Ttllen(Tabcol(Ntab+1))) 250 Format(' Column title is: ',A1/ $ ' Is it ok? (y/n) ',$) Call Instrg('YyNn',4,I) c c If the user doesn't like their choice; c If (I .GT. 2) Then Write(7,300) 300 Format(' Do you wish to stop? (y/n) ',$) Call Instrg('YyNn',4,I) If (I .LE. 2) Then Iflag = 1 Return End If Go To 10 End If c Write(7,350) 350 Format(/' Enter the corresponding column header'/ $ ' (enter a backslash "\"to stop)') End If c c Start reading in the header text c 20 Continue Call Point(Istart,Loop,Tablen) Call Linedt(Table(Istart),Tablen(Loop),80,0) Loop = Loop + 1 Tcolen(Ntab+1) = Tcolen(Ntab+1) + 1 c c Check for a backslash (signifies end of the header) c Do 30 I=1,Tablen(Loop-1) If (Table(I+Istart-1) .EQ. '\' ) Go To 40 30 Continue Go To 20 40 Continue Tablen(Loop-1) = I-1 c c Go here if a backslash is encountered - c see if the user is done c If (I .EQ. 1) Tcolen(Ntab+1) = Tcolen(Ntab+1) - 1 If (Tcolen(Ntab+1) .EQ. 0 .AND. Ntab .GT. 0) Then Iflag = 1 Else If (Ntab .GT. 0) Then c c For data columns; get the number of decimal places c 50 Continue Write(7,450) 450 Format(' Enter the number of decimal places to show: ',$) Call Inchek(0,9,Ndeci(Ntab+1)) c Read(5,*,END=50,ERR=50) Ndeci(Ntab+1) Round(Ntab+1) = 0.0 Icom(Ntab+1) = 0 Isci(Ntab+1) = 0 End If Ntab = Ntab + 1 Iflag = 0 End If c Return End Subroutine Tabspc(Data) c c *********************************************************************** c * * c * Subroutine TABSPC - To get from the user the number of * c * spaces to be put between each of the columns for the table. * c * The current number of spaces are displayed to assist the * c * user in selecting the new number of spaces. * c * * c * Programmed by Steven C. Fisher, DTNSRDC Code 1522 * c * * c *********************************************************************** c Include 'Cmain' Include 'Ctable' c Virtual Data(1) Dimension Isfig(20),Ihedln(20) c c Compute the spacing between the columns c If (Ntab .LT. 3) Return c c Compute the number of spaces between each column c Call Tabwid(Data,Nspac,Ntab,Maxwid) c c Write out the current spacing, and get the new spacing c from the user c Write(7,100) Maxwid 100 Format(/' Currently, the table takes up ',I3,' columns') c 20 Continue c Write(7,200) (I,I+1,Nspac(I),I=1,Ntab-2) 200 Format(' Current number of spaces between each column:'// $ $ (' Between columns',I3,' and',I3,': ',I2,' spaces'/)/ $ ' Do you wish to change the number of spaces? (y/n) ',$) Call Instrg('YyNn',4,J) c If (J .GT. 2) Return c c Getting the new number of spaces c Do 30 I=1,Ntab-2 Write(7,300) I,I+1 300 Format(' Enter the number of spaces between ', $ 'columns',I3,' and',I3,': '$) Call Inchek(1,50,N) If (N .NE. Nspac(I)) Then Maxwid = N - Nspac(I) + Maxwid Write(7,100) Maxwid End If Nspac(I) = N 30 Continue c Go To 20 c End Subroutine SAVTAB(Iflag) c c *********************************************************************** c * * c * Subroutine SAVTAB - To write the current table presentation * c * to a file. * c * * c * Programmed by Steven C. Fisher, DTNSRDC Code 1522 * c * * c *********************************************************************** c Include 'Cmain' Include 'Ctable' Include 'Cflags' c Byte Filnam(25) c c If there is no data, return c c If (Ntab .EQ. 0 ) Return c c Get the file name to store the data c If (Iflag .EQ. 0) Then Write(7,100) 27,27,27 100 Format(1X,1A1,'[H',1A1,'[J'/1X,1A1,'#6 SAVE TABLE PRESENTAT', $ 'ION ROUTINE'//' Enter the file name to store the data '$) Read(5,150,END=98,ERR=98) Filnam 150 Format(25A1) c c Check the file name syntax c Call Chknam(Filnam,Iflg) If (Iflg .EQ. 1) Return c End If c c Compute the filespace needed c Call Point(I,Ntab+1,Tcolen) Call Point(Isave,I,Tablen) Isave = Isave + 100 * Ntab c c Check to see if it is a good file name, and if there is c enough room for the data (Isave is the number of blocks needed c Isave = (Isave + 512)/512 + 2 c If (Iflag .EQ. 0) Call Filchk(Filnam,Isav,Iflg) c c If there is an error with the file name/size/etc, return c If (Iflg .EQ. 1 .AND. Iflag .EQ. 0) Return c If (Iflag .EQ. 0) Then Open(Unit=12,NAME=Filnam,TYPE='NEW',ERR=98,INITIALSIZE=Isave) Else c If (Ntab .LE. 0) Return Open(Unit=12,NAME='PLTPRG.TBL',TYPE='NEW',ERR=98, $ INITIALSIZE=Isave) End If c c Save number of points/column, number of characters/title c c Izero = 0 ! Izero is used to make a space for future improvements c Write(12,225,END=99,ERR=99) Ntab,Iprtln,Izero,Izero 225 Format(16I5) Write(12,225,END=99,ERR=99) (Nspac(I),I=1,Ntab-2) L = 0 Do 20 I=1,Ntab Write(12,230,END=99,ERR=99) Tabcol(I),Tcolen(I),Ndeci(I), $ Icom(I),Isci(I),Round(I),Izero,Izero 230 Format(5I5,G15.7,2I5) Do 10 J=1,Tcolen(I) L = L + 1 Call Point(Istart,L,Tablen) Write(12,250,END=99,ERR=99) (Table(K),K=Istart, $ Istart+Tablen(L)-1) 250 Format(A1) 10 Continue 20 Continue c Close(Unit=12,ERR=99) If (Iflag .NE. 0) Itflag = 0 Return c c Error returns c 98 Continue Iflag = 2 Write(7,300) 300 Format(' WARNING - Error encountered in opening file - ', $ 'routine aborted') Return 99 Continue c Iflag = 2 c Write(7,350) 350 Format(' WARNING - Error encountered in writing data - ', $ 'routine aborted') Close(Unit=12,ERR=101) 101 Continue Return c End