!!$ -*- f90 -*-
! *************************************************************         
!  {\bf giff}                                                           
!                                                                       
! general purpose plot with graphic interactive filtering               
!                                                                       
! orbit9/GNUPLOT version, Pisa, April 1995                              
! *************************************************************         
      PROGRAM giff 
      USE var_precision
      USE fund_const
      USE fir_filters
      USE parameters_giff
      USE prompting
      IMPLICIT REAL(KIND=r_kind)  (a - h, o - z) 
      REAL (KIND=r_kind) xf, yf(2)
! time                                                                  
      DIMENSION tv (nx), x (nx), isel (nx) 
! input data                                                            
      DIMENSION sv (nx, ncox) 
! selected argument                                                     
      DIMENSION crit (nx), y (nx), dy (nx), s (nx) 
! data for filtering, spectra                                           
      DIMENSION amp (nx), tper (nx) 
! control                                                               
      LOGICAL indev, outf 
! system descriptors                                                    
      CHARACTER(3) coox, sysx, unitx, uniy 
      CHARACTER(6) refx 
! names and locations in input file                                     
      CHARACTER(9) num (nplx) 
      DIMENSION iast (nplx) 
! labels                                                                
      CHARACTER ylab * 80, xlab * 80, xlab1 * 80, title * 80 
! filenames                                                             
      CHARACTER spefile * 60, resfile * 60, outfile * 60, postfi * 6 
! *************************************************************   
! device initialisations                                                
      indev = .true. 
      ipag = 0 
      idev=3
      CALL prompt(' device 3=X11 4=Tektr -5=ps file -6=laser pri. -7=color pri.', &
           idev, -7, 6)
! *************************************************************
!  begin again
  999 CONTINUE
! *************************************************************         
! input  mode 
      inter=1                                                          
      CALL prompt(' 1=interactive 0=from giffv.inp', inter, 0, 1) 
      IF (inter.eq.0) then 
         OPEN (4, file = 'giff.inp', status = 'unknown') 
         WRITE ( * , * ) ' reading from file giff.inp' 
         CALL skip (4, 1) 
      ENDIF 
! *************************************************************         
! title                                                                 
      IF (inter.eq.0) then 
         CALL reastr (4, 'title', title) 
         CALL reastr (4, 'postfi', postfi) 
         CALL reaint (4, 'isou', isou) 
      ELSE 
         CALL prompt(' plot title?', title) 
         CALL prompt(' postfix for output files ?', postfi) 
         CALL prompt(' x ax label 1=o8v 2=time 3=each plot 0=no', &
              isou, 0, 3) 
      ENDIF 
! open session log file                                                 
      outfile = 'out.'//postfi 
      CALL remove_spaces (outfile, len) 
      OPEN (10, file = outfile (1:len) , status = 'unknown') 
! x axis labels                                                         
      IF (isou.eq.1) then 
         xlab = 'Time (y) - Orbit8v integration' 
         xlab1 = 'Period (y) - Orbit8v integration' 
      ELSEIF (isou.eq.2) then 
         xlab = 'Time (years)' 
         xlab1 = 'Period (years)' 
      ELSEIF (isou.eq.0) then 
         xlab = ' ' 
         xlab1 = ' ' 
      ENDIF 
! *************************************************************         
!  job input                                                            
      nco = 0 
      nin = 0 
      DO 1 nfile = 1, ncox 
         CALL jobinv (inter, nast, imas, ifo, num, iast, nia, coox,     &
         sysx, refx, unitx, ilce, uniy, ndat)                           
! *************************************************************         
! input data                                                            
         CALL datinp (nast, ifo, iast, nia, ilce, sv, tv, nco, nin,     &
         coox, unitx, uniy, inter, ndat)                                
         CLOSE (1) 
!  input finished?                                                      
         IF (inter.eq.0) then 
            CALL reaint (4, 'imore', imore) 
         ELSE 
            imore=0
            CALL prompt(' more? 1=yes 0=no', imore, 0, 1) 
         ENDIF 
         IF (imore.eq.0) goto 2 
    1 END DO 
! *************************************************************         
! input complete                                                        
    2 WRITE ( *, 150) nin, nco 
      WRITE (10, 150) nin, nco 
  150 FORMAT(' input complete, records=',i6,' with ',i4,' components'/) 
! *************************************************************         
! main loop                                                             
   11 CONTINUE 
! choice of time interval                                               
      CALL timsel (nin, tv, x, isel, nn) 
! *************************************************************         
! choice of argument                                                    
      CALL argsel (nin, nn, nx, nco, x, sv, crit, y, isel, uniy) 
! *************************************************************         
!   choice between plot, filter and spectrum                            
!   menu and function loop                                              
   89 CONTINUE 
      WRITE ( *, 188) 
  188 FORMAT(' 1=plot 2=spectrum 3=filter 4=linfit 5=smooth'/           &
     &   ' 6=out 7=delta,sigma 8=store 9=change device '/               &
     &' 0=change arg -1=stop')                                          
      READ ( *, * ) isp 
      WRITE ( *, 667) isp 
  667 FORMAT(' isp= ',i4) 
      WRITE (10, 667) isp 
! main menu
menu: SELECT CASE(isp)   
      CASE(0) menu  !  change argument
          goto 11 
! **************************************************************        
      CASE(1) menu ! plot of argument    
!  do you want to use the original data, or the transformed ones?  
          call traori (iav, crit, y, s, nn)                            
! device initialization                                                 
         IF (indev) call grafin (idev) 
         indev = .false. 
! new page                                                              
         CALL newpag (idev, ipag) 
! axis labels, if device= printer
         IF (idev.le.0) then 
            IF (isou.eq.3) then  ! if x axis label is not constant
               CALL prompt (' x axis label?', xlab) 
            ENDIF 
! y axis label                                                          
            CALL prompt(' y axis label?', ylab) 
         ENDIF 
! plot
         istyle =1 ! default is continous line
         CALL prompt(' plot style? 0=points 1=line -n=thick points', &
              istyle, -10, 1)
         CALL plotcr (x, s, nn, xlab, ylab, title, idev, istyle) 
         CALL ploflu (idev)
! **************************************************************        
      CASE(2) menu ! spectrum
!  do you want to use the original data, or the transformed ones?       
         call traori (iav, crit, y, s, nn)
         PRINT *, ' period min, max, no. freq. ; 1=print,0=no' 
         READ ( *, * ) pmin, pmax, nfre, iprint 
         WRITE ( *, 668) pmin, pmax, nfre, iprint 
  668 FORMAT(' from ',f12.1,'  to  ',f12.1 ,' no. fr',i5,' ipr',i2) 
         WRITE (10, 668) pmin, pmax, nfre, iprint 
         IF (iprint.ne.0) then 
            spefile = 'spe.'//postfi 
            CALL remove_spaces(spefile, len) 
            OPEN (8, file = spefile (1:len) , status = 'unknown') 
            REWIND 8 
         ENDIF 
         CALL prompt(' 0=ampl  1=sqrt(dens.spett.)', iam, 0, 1) 
         dper = (pmax - pmin) / nfre 
         nfre1 = nfre+1 
         DO  j = 1, nfre1 
            per = pmin + dper * (j - 1) 
            itest = 0 
            CALL peri2 (x, y, nn, per, sp, itest, d0, d1, d2)
            a = d1 * d1 + d2 * d2 
            a = sqrt (a) 
            f = atan2 (d1, d2) 
            cost = + d0 
            IF (iprint.ne.0) write (8, 197) per, sp, a, f, cost 
            tper (j) = per 
            amp (j) = a 
            IF (iam.gt.0) amp (j) = sqrt (sp) 
         END DO 
         IF (iprint.ne.0) close (8) 
! **************************************************************        
! plot of spectrum                                                      
   75    CONTINUE 
! device initialization                                                 
         IF (indev) call grafin (idev) 
         indev = .false. 
! new page                                                              
         CALL newpag (idev, ipag) 
         IF (idev.lt.0) then 
! x axis label                                                          
            IF (isou.eq.3) then 
               CALL prompt (' x axis label?', xlab1) 
            ENDIF 
! y axis label                                                          
            CALL prompt( ' y axis label?', ylab) 
         ENDIF 
! plot                                                                  
         istyle = 1 
         CALL plotcr (tper, amp, nfre1, xlab1, ylab, title, idev, istyle)
         CALL ploflu (idev) 
!  what next?                                                           
         IF (iprint.ne.0) close (8) 
            ndev=0 !default
            CALL prompt(' change device? -5=ps file -6=laser 3=X11 4=TeK 0=as before', &
                 ndev, -7, 6)
         IF (ndev.ne.0) then 
            CALL grafcl (idev) 
            indev = .true. 
            ipag = 0 
            idev = ndev 
            GOTO 75 
         ENDIF 
! **************************************************************
      CASE(3) menu ! filtering : Ferraz-Mello single frequency filter 
!  do you want to use the original data, or the transformed ones?       
         call traori (iav, crit, y, s, nn)   
         WRITE ( *, 194) ! choice of frequency to be filtered out      
  194 FORMAT   (' output: period,rel.sp.dens.,amp,phas,const/sig'/) 
         WRITE (10, 194) 
   70    CONTINUE 
         WRITE ( *, 193) 
  193 FORMAT   (' period, itest(1=filter,0=no); per<0 to stop' ) 
         READ ( *, * ) per, itest 
         IF (per.le.0.d0) goto 89 
         CALL peri2 (x, s,  nn, per, sp, itest, d0, d1, d2) 
         a = d1 * d1 + d2 * d2 
         a = sqrt (a) 
         f = atan2 (d1, d2) 
         cost = d0 
         WRITE ( *, 197) per, sp * 100, a, f, cost, itest 
  197 FORMAT   (f14.4,f8.4,1p,d13.6,0p,f9.5,1p,d13.6,i2) 
         WRITE (10, 197) per, sp, a, f, cost, itest 
         IF (itest.ne.0) then 
            sm=sum(s(1:nn))/nn
            sig = 0.d0 
            DO  i = 1, nn 
               sig = sig + (s (i) - sm) **2 
               y (i) = s (i) 
            END DO 
            sig = sqrt (sig / nn) 
            WRITE ( *, 198) sig 
  198 FORMAT      (1p,d15.5) 
            WRITE (10, 198) sig 
         ENDIF 
         GOTO 70 
! *************************************************************         
      CASE(4) menu !  linear fit
!  do you want to use the original data, or the transformed ones?       
         call traori (iav, crit, y, s, nn)
         CALL linfi3 (x, s, rate, rms, cost, dy, nn) 
         WRITE ( *, 941) rate, cost, rms, iav 
  941 FORMAT   (' rate=',1p,d18.10,'  cost=',d15.7,'  rms=',d12.4,i5) 
         WRITE (10, 941) rate, cost, rms, iav 
         y (1:nn) = dy (1:nn)
! *************************************************************         
      CASE(5) menu !  digital smoothing
!  do you want to use the original data, or the transformed ones?       
         call traori (iav, crit, y, s, nn)
         WRITE ( * , * ) ' sampling ratio?' 
         READ ( *, * ) isamp 
         WRITE (10, 959) isamp 
  959 FORMAT   (' digital smoothing, decimation ',i4) 
!  coefficienti del filtro  
         CALL filope(isamp,1)
!  passo assunto costante
         fstep= x(nin/2+1) - x(nin/2)
! filter half length
         hlen=fstep*(nfil-1)/2 
!  digital low frequency pass filter                                    
         nskip = 0 
         j = 0 
         DO 950 i = 1, nn 
            ndat=1
            CALL filter (x(i), s(i:i), ndat, nskip, fstep, xf, yf(1:1), outf)
            IF (outf) then 
               j = j + 1 
               y (j) = yf(1)
               x (j) = xf
            ENDIF 
  950    END DO 
         nn = j 
         WRITE ( *, 953) nn 
  953 FORMAT   (' available data points no.', i5/                       &
     &    '  warning: times are changed; do not use orig.data')  
! *************************************************************         
      CASE(6) menu !  output data as they are now in file res.postfi  
!  do you want to use the original data, or the transformed ones?       
         call traori (iav, crit, y, s, nn)
         resfile = 'res.'//postfi 
         CALL remove_spaces (resfile, len) 
         OPEN (11, file = resfile (1:len) , status = 'unknown') 
         DO 963 i = 1, nn 
          WRITE (11, 962) x (i), s (i)
  963    END DO           ! ADDED
  962 FORMAT   (f15.3,e18.9) 
         CLOSE (11) 
! *************************************************************  
      CASE(7) menu !  average and sigma ,max and min 
!  do you want to use the original data, or the transformed ones?       
         call traori (iav, crit, y, s, nn)
         swy = 0.d0 
         smax = maxval(s(1:nn))
         smin = minval(s(1:nn))
         swy = sum(s(1:nn))
         ym = swy / nn 
         del = smax - smin 
         WRITE ( *, 181) ym, smax, smin, del 
  181 FORMAT(' av= ',1p,d15.7,' max=',d15.7,' min=',d15.7,' del=',d15.7)
         WRITE (10, 181) ym 
         sig = 0.d0 
         ii = 0 
         DO  i = 1, nn 
            sig = sig + (s (i) - ym) **2 
         END DO 
         sig = sqrt (sig / nn) 
         WRITE ( *, 182) sig 
  182 FORMAT   (' sigma=',1p,d13.5) 
         WRITE (10, 182) sig 
! *************************************************************  
      CASE(8) menu !  store current argument                                   
         IF (nco.eq.ncox) then 
      WRITE ( * ,  * ) ' not enough room, already ', ncox, ' components' 
            GOTO 89 
         ELSEIF (nin.ne.nn) then 
            WRITE ( * , * ) ' this is not allowed, you can only store' 
            WRITE ( * , * ) ' over the entire time span' 
            GOTO 89 
         ENDIF 
!  do you want to use the original data, or the transformed ones?       
         call traori (iav, crit, y, s, nn)
         nco = nco + 1 
         sv (1:nn, nco) = s (1:nn) 
         WRITE ( * , * ) ' stored in column ', nco 
! *************************************************************         
      CASE(9) menu !  change device                                            
         CALL prompt(' change device? -5=ps file -6=laser 3=X11 4=Tek 0=as before', &
              idev, -7, 6)
         IF(idev.ne.0)THEN
            CALL grafcl (idev) 
            indev = .true.   
            ipag = 0   
         ENDIF        
! **************************************************************        
      CASE(-1) menu ! close and/or change file
         CALL PROMPT(' change input file? 1=yes 0=stop', &
              ibegin, 0, 1)
         CLOSE(10)
         IF(ibegin.eq.1)THEN
            goto 999
         ENDIF 
         STOP 
! *************************************************************         
! input error                                                           
      CASE DEFAULT  menu
         WRITE ( * , * ) ' option ', isp, ' unknown' 
      END SELECT menu
      GOTO 89 
      END PROGRAM giff                              
! *************************************************************         
!   {\bf argsel} choice of argument                                     
      SUBROUTINE argsel (nin, nn, nx, nco, x, sv, crit, y, isel, uniy)
      USE var_precision
      USE fund_const
      USE prompting
      IMPLICIT real(kind=r_kind) (a - h, o - z)
      PARAMETER (ncox = 100) 
      CHARACTER(3) uniy 
      LOGICAL error
      DIMENSION x (nn), crit (nn), sv (nx, nco), isel (nin), y (nn),    &
      k (ncox)
   40 CONTINUE  ! just to pass control by later GOTO's
      CALL prompt(' which data to plot? 0=comb -1=two arg -2=ratio & 
         & -3 polar radius -4 polar angle -5=cart from polar', js, -5, nco)
      WRITE ( *, 600) js 
  600 FORMAT(' data choice: ',i5) 
      WRITE (10, 600) 
! ************************************
      IF(js.lt.0) THEN  ! plot of two coordinates
!  to handle a couple of coordinates (instead of time and one coordinate
!                                                                       
!  select the two coordinates 
         WRITE ( * , * ) ' two indices?' 
         j1=0 ; j2=0 !no defaults, the range may have changed
         CALL prompt(' first coordinate ', j1, 1, nco)
         CALL prompt(' second coordinate ', j2, 1, nco)
axes:    SELECT CASE(js) 
! ************************************      
         CASE (-1) axes ! plot two coordinates as they are                    
               WRITE (10, 1410) j1, j2 
 1410          FORMAT(' 2 components, no.',2i5) 
               ii = 0 
               DO i = 1, nin 
                  IF (isel (i) .eq.1) then 
                     ii = ii + 1 
                     x (ii) = sv (i, j1) 
                     crit (ii) = sv (i, j2) 
                  ENDIF 
               END DO 
! ************************************ 
         CASE(-2) axes  !  ratio of two coordinates    
            WRITE (10, 1411) j1, j2 
 1411       FORMAT (' ratio of ',i5,'  and ',i5) 
            ii = 0 
            izer = 0 
            DO  i = 1, nin 
               IF (isel (i) .eq.1) then 
                  ii = ii + 1 
                  IF (sv (i, j2) .eq.0.d0) then 
                     izer = izer + 1 
                     ii = ii - 1 
                  ELSE 
                     crit (ii) = sv (i, j1) / sv (i, j2) 
                  ENDIF 
               ENDIF 
            END DO 
            IF (izer.gt.0) write ( * , * ) ' zero divide no. ', izer 
! ************************************
         CASE( -3) axes  !  polar radius 
            WRITE (10, 1412) j1, j2 
 1412       FORMAT (' polar radius of ',i5,'  and ',i5) 
            ii=0
            DO  i = 1, nin 
               IF (isel (i) .eq.1) then 
                  ii = ii + 1 
                  crit (ii) = sqrt (sv (i, j2) **2 + sv (i, j1) **2) 
               ENDIF 
            END DO 
! ************************************
         CASE( -4) axes  !  polar angle, in the current unit for angles 
            WRITE (10, 1413) j1, j2, uniy 
 1413 FORMAT       (' polar angle of ',i5,'  and ',i5, 'unit= ',a3) 
            CALL convun_giff ('RAD', uniy, uconv, urev, error) 
! angle update - sostituire con type angle
            ng = 0 
            crv = 0.d0 
            ns = 0 
            ii=0
            DO 394 i = 1, nin 
               IF (isel (i) .eq.1) then 
                  ns = ns + 1 
                  ii = ii + 1 
                  crr = atan2 (sv (i, j2), sv (i, j1) ) 
                  IF (crr.lt.crv - pig) then 
                     ng = ng + 1 
                  ELSEIF (crr.gt.crv + pig) then 
                     ng = ng - 1 
                  ENDIF 
                  crit (ii) = crr + pig * ng * 2.d0 
                  crv = crr 
               ENDIF 
  394       END DO 
            crit (1:ns) = crit (1:ns) * uconv 
! ************************************
         CASE(-5) axes ! cartesian coordinates from polar ones
            WRITE (10, 1414) j1, j2, uniy 
 1414 FORMAT       (' cartesian coord of ',i5,' as radius '/            &
     &         10x,i5, ' as angle, unit= ',a3)                          
            CALL convun_giff ('RAD', uniy, uconv, urev, error) 
            ns = 0
            ii=0 
            DO  i = 1, nin 
               IF (isel (i) .eq.1) then 
                  ns = ns + 1 
                  ii = ii + 1 
                  crr = sv (i, j2) / uconv 
                  rr = sv (i, j1) 
                  crit (ii) = rr * sin (crr) 
                  x (ii) = rr * cos (crr) 
               ENDIF 
            END DO 
         CASE DEFAULT axes ! unknown negative js
            PRINT *, ' Unknown option '
            GOTO 40
         END SELECT axes
! ***********************************************                             
      ELSEIF (js.eq.0) then 
!  linear combination, with integer arguments, of the coordinates 
         WRITE ( *, 141) nco 
  141 FORMAT   (i4,' integers?') 
         DO  j = 1, nco 
            READ ( *, * ) k (j) 
         END DO
         WRITE ( *, 142) (k (i), i = 1, nco) 
         WRITE (10, 142) (k (i), i = 1, nco) 
  142    FORMAT(' combination,  coeff. '/(10i4)) 
         ii = 0 
         DO  i = 1, nin 
            IF (isel (i) .eq.1) then 
               ii = ii + 1 
               crit(ii)=dot_product(k(1:nco),sv(i,1:nco))
            ENDIF 
         END DO 
! ***********************************************
      ELSEIF (js.le.nco) then 
!  plot time and one coordinate    
         WRITE (10, 1440) js 
 1440    FORMAT(' component no. ',i4) 
         ii = 0 
         DO  i = 1, nin 
            IF (isel (i) .eq.1) then 
               ii = ii + 1 
               crit (ii) = sv (i, js) 
            ENDIF 
         END DO 
      ELSE 
         WRITE ( *, 143) j, nco 
  143    FORMAT(' component ',i4,' not available, maximum ',i4) 
         GOTO 40 
      ENDIF 
! *************************************************************         
! argument has been selected                                            
!                                                                       
!  average and sigma                                                    
      ym = sum(crit(1:nn))/nn
      WRITE ( *, 181) ym 
  181 FORMAT( ' average= ',1p,d20.11) 
      WRITE (10, 181) ym 
      y(1:nn)=crit(1:nn)-ym
      sig=sqrt(dot_product(y(1:nn),y(1:nn))/nn)
      WRITE ( *, 182) sig 
  182 FORMAT(' sigma=',1p,d13.5) 
      WRITE (10, 182) sig 
      RETURN 
      END SUBROUTINE argsel                         
! *************************************************************         
!  {\bf timsel} time interval selection                                 
      SUBROUTINE timsel (nin, tv, x, isel, nn) 
      USE var_precision
      USE prompting
      IMPLICIT real(kind=r_kind) (a - h, o - z)      
! non occorre      DOUBLE PRECISION dta,dtz
      DIMENSION tv (nin), x (nin), isel (nin) 
! chose criterium                                                       
   50 WRITE ( *, 151) nin, tv (1), tv (nin) 
  151 FORMAT(i6,' data points, from t=',f16.6,'  to t=',f16.6)
      iin=0
      CALL prompt (' time int. 0=all 1=by no. point 2=by time ?', &
           iin, 0, 2)                 
      IF (iin.eq.0) then 
         ia = 1 
         iz = nin 
         nn = nin 
         isel (1:nn) = 1 
      ELSEIF (iin.eq.1) then 
         WRITE ( * , * ) ' initial an final data point no. ?' 
         ia=1 !default
         CALL prompt(' initial ',ia, 0, nin)
         iz=nin !default
         CALL prompt(' final ',iz, ia, nin)
         nn = iz - ia + 1  
         isel (1:nin) = 0  
         isel (ia:iz) = 1 
      ELSEIF (iin.eq.2) then 
         tmin=minval(tv(1:nin))
         tmax=maxval(tv(1:nin))
         WRITE ( * , * ) ' initial and final time?'
         ta=tmin !default 
         CALL prompt(' initial ',ta, tmin, tmax)
         tz=tmax !default
         CALL prompt(' final ',tz, ta, tmax)
         nn = 0 
         DO 51 n = 1, nin 
            IF (tv (n) .le.tz.and.tv (n) .ge.ta) then 
               nn = nn + 1 
               isel (n) = 1 
            ELSE 
               isel (n) = 0 
            ENDIF 
   51    END DO 
      ELSE 
         WRITE ( *, 154) iin 
  154 FORMAT   (i4,' option not supported') 
         GOTO 50 
      ENDIF 
      IF (nn.gt.0) then 
         ns = 0 
         DO 522 i = 1, nin 
            IF (isel (i) .eq.1) then 
               ns = ns + 1 
               x (ns) = tv (i) 
            ENDIF 
  522    END DO 
      ELSE 
         WRITE ( * , * ) ' no interval selected' 
         GOTO 50 
      ENDIF 
!                                                                       
      WRITE (10, 1520) nn, x (1), x (nn) 
 1520 FORMAT(i6,' data points, time from ',f15.5,'  to  ',f15.5) 
      RETURN 
      END SUBROUTINE timsel                         
! ************************************************************          
!  {\bf traori} choice between original and transformed data            
      SUBROUTINE traori (iav, crit, y, s, nn) 
      USE var_precision
      USE prompting
      IMPLICIT real(kind=r_kind) (a - h, o - z) 
      DIMENSION crit (nn), y (nn), s (nn) 
      iav=0 !default: original data
      CALL prompt(' 1=transformed data 0=original ?', iav, 0, 1) 
      WRITE (10, 601) iav 
  601 FORMAT(' iav=',i4) 
      IF (iav.eq.0) then 
        s (1:nn) = crit (1:nn) 
      ELSEIF (iav.eq.1) then 
        s (1:nn) = y (1:nn) 
      ENDIF 
      RETURN 
      END SUBROUTINE traori                         


