Subroutine QUIKPL(Data) c c *********************************************************************** c * * c * Subroutine QUIKPL - To create a "Quick" plot presentation * c * to plot out the users data. This routine asks a minimum of * c * information from the user - the number of curves to plot, * c * the X and Y column pairs, and the corresponding curve fit. * c * The rest of the parameters usually needed to create a * c * presentation are set by QUIKPL. If a user currently has * c * a plot presentation, it is not destroyed by this routine. * c * The user created presentation is stored as PLTPRG.PLT. * c * (PLTPRG always reads PLTPRG.PLT along with PLTPRG.DAT and * c * PLTPRG.TBL when it restarts). * c * * c * Programmed by Steven C. Fisher, DTNSRDC Code 1522 * c * * c *********************************************************************** c Include 'Cmain' Include 'Cflags' c Virtual Data(1) Dimension Ixcrv(20),Iycrv(20),Ichos(20),Njunct(20),Xloc(60) Character Com*40,File_Name*60 Integer Choice Integer*4 Isave Character*7 Xax,Yax Character*22 Label c Data Xax,Yax/' X Axis',' Y Axis'/ Data R9,R6p5,I4,Rp01,Rp38,R7p5,R1p1,Zero,Rone/9,6.5,4,0.01,0.38, $ 7.5,1.1,0.0,1.0/ Data Izero,Ione/0,1/ c c See if there is any data to interpolate c If (Ncols .EQ. 0) Then Write(6,110) 110 Format(' ERROR - No Data in the program - routine aborted') Return Endif c c Isave contains the number of bytes needed in file c PLTPRG.PLT to hold the plot presentation c Isave = 901 c c Input the curve to interpolate c Write(6,100) 27,27,27 100 Format(1X,1A1,'[H',1A1,'[J'/1X,1A1,'#6 DATA PLOTTING ', $ 'ROUTINE'/) Write(6,150) 150 Format(' Enter the number of curves to plot') Call Indata(1,'I',Ncurv) If (Ncurv .LE. 0) Return c Maxpt = 0 Do 25 II = 1,Ncurv Label = ' ' 15 Continue Write(6,250) II,Label 250 Format(' Enter the column numbers corresponding to'/ $ ' the X and Y values for curve no. ',I3,3X,A) Call Inpt(X,N,Iflag,2) If (Iflag .EQ. 1) Return If (N .EQ. 0) Go To 15 Ix = X(1) If (N .EQ. 1) Then Call Inpt(X,N,Iflag,2) If (Iflag .EQ. 1) Return If (N .EQ. 0) Go To 15 Iy = X(1) Else Iy = X(2) End If 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,6) If (Choice .EQ. 1) Then Label = '(or "S" to stop)' Go To 15 End If c c See if only one point is in each column (error condition) c If (Npoint(Ix) .EQ. 1) Then Write(6,300) 300 Format(/' ERROR - Only one point in each column'/) Label = '(or "S" to stop)' Go To 15 End If c c See if the user likes his column choices c Write(6,350) 350 Format(' Are these column numbers correct? (y/n) '$) Call Instrg('YyNn',4,Choice) If (Choice .GT. 2) Then Label = '(or "S" to stop)' Go To 15 End If Ixcrv(II) = Ix Iycrv(II) = Iy c c Get the data curve fit type c Write(6,400) 400 Format(' Enter the curve fit to interpolate:'// $ ' 0) Draw points only'/ $ ' 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 fit - type A'/ $ ' 8) Least squares cubic spline fit - type B'// $ ' Enter 0 - 8 '$) Call Inchek(0,8,Ichos(II)) c c If a Least Squares is selected, get the fit order: c If (Ichos(II) .EQ. 2) Then Write(6,425) 425 Format(' Enter the order of the polynomial fit (Cubic=3)', $ ' ',$) Call Inchek(1,20,J) Ichos(II) = -J End If c If (Ichos(II) .GE. 7) Then c c Least squares cubic spline fit c Label = ' ' Nl = 1 40 Continue Write(6,430) 58+Ichos(II),Label 430 Format(/' Least Sqaures Cubic Spline Fit Type ',1A1,':'/ $ ' Enter the column number containing the least'/ $ ' squares spline X junctures ',A,' '$) Call Inpt(X,N,Iflag,2) If (Iflag .EQ. 1) Then Label = ' ' Go To 15 End If If (N .LT. 1) Go To 40 If (X(1) .LT. 1 .OR. X(1) .GT. Ncols) Then Write(6,377) Ncols,7 377 Format(' ERROR - Column number must be between 1', $ ' and',I4,1A1) Label = '(or "S" to redo curve)' Nl = 22 Go To 40 End If Icol = X(1) Call Name(Titl,Icol,Title,Ttllen) Write(6,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) Then Label = '(or "S" to redo curve)' Nl = 22 Go To 40 End If c c See if we have too many points in this junction column c If (Npoint(Icol) .GT. 60) Then Write(6,385) 7 385 Format(' ERROR - Too many points in this column - ', $ 'Try again',1A1) Label = '(or "S" to redo curve)' Nl = 22 Go To 40 End If c Call Getvec(Data,Npoint,Icol,Xloc,Nxloc) 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(6,395) Xloc(I),Xloc(I+1) 395 Format(' ERROR - There are no datapoints between the'/ $ ' junction points X= ',G12.5, $ ' and X= ',G12.5/) Label = '(or "S" to redo curve)' Nl = 22 Go To 40 End If 60 Continue Njunct(II) = Icol End If c c Compute the maximum and minimum values for X and Y c Call Getvec(Data,Npoint,Ix,X,Npt) Call Getvec(Data,Npoint,Iy,Y,Npt) If (II .EQ. 1) Then Xmin = X(1) Xmax = X(1) Ymin = Y(1) Ymax = Y(1) End If c Isave = Isave + Npt * 26 + 40 c If (Ichos(II) .GE. 7) Isave = Isave + Npoint(Njunct(II))*13+15 c c To find the maximum and minimum Y values for the plot, we c play a little game. We look at the size of the difference c between the Y values near the maximum and minimum Y values. c The reason for this is that when we fit a curve to the data, c the curve may go higher (or lower) than the maximum (or minimum) c Y data value. So, to ensure that these curve peaks and troughs c will not be cut off on the plot, and since the size of the peak c of the curve relative to the nearest point is related to the c difference in the neighboring Y values, we did this: c Do 23 I=1,Npt c If (I .EQ. 1) Then Delta = Abs((Y(2) - Y(1)) * 0.4) Else If (I .LT. Npt) Then D1 = Abs(Y(I) - Y(I-1)) D2 = Abs(Y(I+1) - Y(I)) If ( D1 .GT. D2 * 3.0) D2 = D1 If ( D2 .GT. D1 * 3.0) D1 = D2 Delta = (D1 + D2) * 0.2 Else Delta = Abs((Y(Npt) - Y(Npt-1)) * 0.4) End If c If (Xmin .GT. X(I)) Xmin = X(I) If (Xmax .LT. X(I)) Xmax = X(I) If (Ymin .GT. Y(I) - Delta) Ymin = Y(I) - Delta If (Ymax .LT. Y(I) + Delta) Ymax = Y(I) + Delta 23 Continue Maxpt = Max0(Maxpt,Npt) 25 Continue Xmin = Xmin - 0.02 * (Xmax - Xmin) Xmax = Xmax + 0.02 * (Xmax - Xmin) c c Get the increments for X and Y axis: c Call Getax(Xmin,Xmax,Xmn,Xmx,Xmajor,Xminor) Call Getax(Ymin,Ymax,Ymn,Ymx,Ymajor,Yminor) c c If the number of points > 25, see if we plot symbols c If (Maxpt .GT. 25) Then Write(6,500) 500 Format(' Do you wish to plot symbols at the datapoints?', $ ' (y/n) ',$) Call Instrg('YyNn',4,Isymb) Isymb = 2-((Isymb+1)/2) Else Isymb = 1 End If c c Get the device to plot to c Call Pltdev(Idev,1,Ibatch) If (Ibatch .GT. 0) Then Iax = 1 Else Iax = 0 End If c c Isave = Isave/512 + 4 c c Get the scratch file name and detach command c Call Getusr(File_Name,Com) c c Open the file to store the data c Open(Unit=12,Name=File_Name,Type='NEW',ERR=97) Write(12,265,ERR=95) Idev c c Xlen,Ylen,Textsz,Axchst,Axthik,Spchar c Write(12,230,ERR=95) R9,R6p5,I4,Iax,Rp01,Ione c c Xpmin,Xpmax,Ypmin,Ypmax c Write(12,270,ERR=95) 0.0,9.0,-1.0,6.7 230 Format(2F5.2,2I4,F5.2,I4,2F5.2) c c Xnumbr,Ynumbr,Axtypx,Axtypy,Ticdir,Atxclr c Write(12,265,ERR=95) Ione,Ione,Ione,Ione,Ione,Ione c c Xmin,Xmax,Xinc c Write(12,270,ERR=95) Xmn,Zero,Xmx,Zero,Xminor,Zero c c Ymin,Ymax,Yinc c Write(12,270,ERR=95) Ymn,Zero,Ymx,Zero,Yminor,Zero c c Grid, # labels & Legends c c Xgrdmj,Xgrdmr,Ygrdmj,Ygrdmr c Write(12,270,ERR=95) Zero,Zero,Zero,Zero c c Thikmj,Thikmr,Grdcmj,Grdcmr c Write(12,245,ERR=95) Zero,Zero,Izero,Izero 245 Format(2F5.1,2I4) c c Quick plot index (=1),Axcolr,Lecolr,Lacolr c Write(12,265,ERR=95) Ione,Ione,Ione,Ione Write(12,270,ERR=95) Zero,Zero,Zero c c Axis labels c Write(12,255,ERR=95) Xax,Yax 255 Format(A7//A7/) c c NLab,Lachst,Ntxbl,Nleg,Lechst,Legtit,Labcrd c Write(12,265,ERR=95) 0,0,0,0,0,0,0,0 c c Ncurve,Cvthik,Symbsz c Write(12,*,ERR=95) Ncurv,Rp01,Rone,Izero 260 Format(I4,2F5.3,I4) Do 70 I=1,Ncurv c J = I If (J .GT. 10) J = J - 10 Isymbl = 1 If (Isymb .EQ. 0 .AND. Ichos(I) .NE. 0) Isymbl = 0 Write(12,265,ERR=95) Npoint(Ixcrv(I)),J,Isymbl*I, $ Ichos(I),Ione,Ione,Ione,Ixcrv(I),Iycrv(I),Izero 265 Format(10I5) 70 Continue c Do 80 I=1,Ncurv c Call Point(Ipoint,Ixcrv(I),Npoint) Write(12,270,ERR=95) (Data(J),J=Ipoint, $ Npoint(Ixcrv(I)) + Ipoint-1) 270 Format(6G13.5) Call Point(Ipoint,Iycrv(I),Npoint) Write(12,270,ERR=95) (Data(J),J=Ipoint, $ Npoint(Iycrv(I)) + Ipoint-1) c If (Ichos(I) .GE. 7) Then c c Write in the control point info for the LSQ spline c Write(12,265,ERR=95) Npoint(Njunct(I)),Njunct(I) Call Point(Ipoint,Njunct(I),Npoint) Write(12,270,ERR=95) (Data(J),J=Ipoint, $ Npoint(Njunct(I)) + Ipoint-1) End If 80 Continue Close(Unit=12) c c See if we can send it as a detached job for the HP plotter c If (Ibatch .EQ. 2) Then Call STDET(Idev,Ierr) If (Ierr .EQ. 0) Then Write(6,600) 600 Format(' Plot is being done as a batch job') Return End If Ibatch = 1 End If c c Store the current presentation on disk c If (Ipflag .EQ. 1) Then Open (Unit=12,NAME='PLTPRG.PLT',TYPE='NEW',Initialsize=7, $ ERR=97) c Write(12,265,ERR=95) Idev Call Wpldat(Data,Iflag,0) If (Iflag .EQ. 1) Return Ipflag = 0 End If c c Store the current data on disk c If (Idflag .EQ. 1) Then Call Wrtnow(Data,Iflag) If (Iflag .EQ. 1) Go To 97 Idflag = 0 End If c c Set flag indicating that you have saved the data c c If (Itflag .EQ. 1) Then Iflag = 1 Call Savtab(Iflag) If (Iflag .EQ. 2) Go To 97 Itflag = 0 End If c c Jump to the plotting program c If (Ibatch .GT. 0) Write(6,650) 650 Format(' Computing plot. Please wait...'/) Call Runcmd(Com) Call Exit c 95 Continue Write(6,700) 700 Format(' File error encountered - routine aborted') Return c 97 Continue Write(6,800) 800 Format(/' ERROR - Insufficient file space for scratch file - ', $ 'routine aborted'/) Return End