Subroutine Initrp(Data,Ix,Iy,Npt,Xvals,Numbr,Xloc,Nxloc, $ Icheck,Iorder,Iflag) c c *********************************************************************** c * * c * Subroutine Initrp - To enable the user to interpolate * c * values from data currently inside PLTPRG. The user can * c * use straight line, cubic spline, parametric spline * c * (3 types), or a least-squares-fit polynomial to fit the * c * existing data, for interpolation and extrapolation. The * c * interpolated values can be stored back as data, if desired. * c * * c * Programmed by Steven C. Fisher, DTNSRDC Code 1522 * c * * c *********************************************************************** c Include 'Cmain' c Virtual Data(1) Dimension Xvals(1),Xloc(1) Integer Choice c c See if there is any data to interpolate c Iflag = 0 c If (Ncols .EQ. 0) Then Write(7,110) 110 Format(' ERROR - No Data in the program - routine aborted') Iflag = 1 Return End If c 10 Continue c c Input the curve to interpolate c Write(7,100) 27,27,27 100 Format(1X,1A1,'[H',1A1,'[J'/1X,1A1,'#6 DATA INTERPOLATION ', $ 'ROUTINE'/) 15 Continue Write(7,105) 105 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,Choice,7) If (Choice .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,120) 120 Format(' Only one point in each column - routine aborted') Iflag = 1 Return End If c c See if the user likes his column choices c Write(7,150) 150 Format(' Are these column numbers correct? (y/n) '$) Call Instrg('YyNn',4,Choice) If (Choice .GT. 2) Then Write(7,155) 155 Format(' Do you wish to stop? (y/n) '$) Call Instrg('YyNn',4,Choice) If (Choice .GT. 2) Go To 10 Iflag = 1 Return End If c c Determine the type of values to interpolate will be input c Write(7,250) 250 Format(' Enter Either'// $ ' 1) Specific X values'/ $ ' 2) Minimum, maximum, and increment values'/ $ ' 3) X values stored in a column'// $ ' Enter 1-3 '$) Call Inchek(1,3,Ich) c c Enter the values to interpolate c If (Ich .EQ. 1) Then c c Specific X values c Write(7,300) 300 Format(' Enter the X values - enter "S" to stop') Call INPT(Xvals,Numbr,If,1) c Else If (Ich .EQ. 2) Then c c Minumum, maximum, and increment c Write(7,350) 350 Format(' Enter the minimum, maximum, and increment X values') Call Indata(3,'RRR',Xmin,Xmax,Xinc) c c Compute the X values c Numbr = 1 + Abs((Xmax - Xmin) / Xinc) If (Numbr .GT. Nptmax) Numbr = Nptmax Do 25 I=1,Numbr Xvals(I) = Amin1(Xmin,Xmax) + Abs((I-1) * Xinc) 25 Continue Else c c X values stored in a column c 30 Continue Write(7,375) 375 Format(' Enter the column number containing the X values to', $ ' interpolate '$) Call Inchek(1,Ncols,Icol) Call Name(Titl,Icol,Title,Ttllen) Write(7,380) (Titl(I),I=1,Ttllen(Icol)) 380 Format(' Column title is:'/1X,A1/ $ ' Is it correct? (y/n) ',$) Call Instrg('YyNn',4,I) If (I .GT. 2) Go To 30 Call Getvec(Data,Npoint,Icol,Xvals,Numbr) c 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) Then Write(7,390) 390 Format(' Sort the X values to interpolate at? (y/n) ',$) Call Instrg('YyNn',4,I) If (I .LT. 3) Call Sort(Xvals,Y,Numbr) End If c c Get the data curve fit type c Write(7,400) 400 Format(' Enter the curve fit to interpolate:'// $ ' 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 segments'/ $ ' 8) Least squares cubic spline fit (type A)'/ $ ' 9) Least squares cubic spline fit (type B)'// $ ' Enter 1 - 9 '$) Call Inchek(1,9,Icheck) c c Get the X and Y values stored in the columns c Call Getvec(Data,Npoint,Ix,X,Npt) c c Sort the data (unless a parametric spline of type c c is selected) c If (Icheck .NE. 6) Then Call Getvec(Data,Npoint,Iy,Y,Npt) Call Sort(X,Y,Npt) Else Call Getvec(Data,Npoint,Iy,Y,Npt) End If c If (Icheck .EQ. 2) Then c c Least squares fit c Write(7,450) 450 Format(' Enter the degree of the polynomial (3 = cubic) '$) Call Inchek(1,20,Iorder) End If c If (Icheck .EQ. 7) Then Write(7,450) Nmx = Min0(16,Npt-1) Call Inchek(1,Nmx,Iorder) Write(7,500) 500 Format(' Enter the number of points fit by each '/ $ ' least squares polynomial segment ',$) Nmx = Min0(Npt,Nptmax) Nm = Max0(3,Iorder+1) Call Inchek(Nm,Nmx,N) Write(7,550) 550 Format(' Do you want the points weighted (y/n) ',$) Call Instrg('YyNn',4,Ichose) If (Ichose .LT. 3) Then Write(7,600) 600 Format(' Select the weight: '// $ ' 1) No additional weight (weight of 1)'/ $ ' 2) Weight of 2'/ $ ' 3) Weight of 5'/ $ ' 4) Weight of 25'/ $ ' 5) Weight of 100'/ $ ' Enter 1-5 ',$) Call Inchek(1,5,Iwt) Else Iwt = 1 End If Iorder = Iorder*1000 + N*10 + (Iwt-1) End If c If (Icheck .GE. 8) Then c c Least squares cubic spline fit c 35 Continue Write(7,650) 650 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,300) Call INPT(Xloc,Nxloc,If,1) Else If (Ich .EQ. 2) Then 40 Continue Write(7,375) Call Inchek(1,Ncols,Icol) Call Name(Titl,Icol,Title,Ttllen) Write(7,380) (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,685) 685 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,700) Xloc(I),Xloc(I+1) 700 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 Return End