      PROGRAM specwork
      IMPLICIT NONE
      INTEGER nspec,column
      INTEGER i,j,n,numdata,nline,step,pp
      INTEGER first(150),last(150),nc,ps
      REAL wl(15),time(150),c(10)
      REAL vmin,vmax,grmin,grmax,t0
      CHARACTER opt
      CHARACTER*64 spec(15,150)
      CHARACTER*16 title,wlname(15),pgname(15),vps(3),vpp(2)
      CHARACTER*16 tblfile,datfile,timefile
      CHARACTER*16 plotmode
      PARAMETER(tblfile='specwork.tbl')
      PARAMETER(datfile='specwork.dat')
      PARAMETER(timefile='specwork.time')
      COMMON /param/ title,nline,nspec,wlname,wl,pgname,time
      plotmode='/XW'
      vmin=-400.
      vmax=400.
      vps(1)='Line Contour'
      vps(2)='Shading'
      vps(3)='Line & Shading'
      vpp(1)='Main'
      vpp(2)='Correlation'
      ps=1
      pp=1
      nc=5
      grmin=0.6
      grmax=1.0
      WRITE(*,*)'Specwork - Astronomical Spectral Analysis (96/06/19)'
      WRITE(*,*)'Stockholm Observatory'
      WRITE(*,*)'Written by'
      WRITE(*,*)'Sten Andersson sten@astro.su.se'
      WRITE(*,*)'Erik Gullbring erik@astro.su.se'
      WRITE(*,*)'Reading ',datfile
      OPEN(UNIT=50,FILE=datfile)
      CALL ReadParam
      CLOSE(50)
      WRITE(*,*)'Reading ',tblfile
      OPEN(UNIT=50,file=tblfile)
      DO i=1,nline
         DO j=1,nspec
            READ(50,5000)spec(i,j)
         ENDDO
      ENDDO
      CLOSE(50)
      WRITE(*,*)'Reading ',timefile
      OPEN(UNIT=50,file=timefile)
      DO i=1,nspec
         READ(50,*)time(i)
         IF (i.EQ.1)t0=time(i)
         time(i)=(time(i)-t0)*24
      ENDDO
      CALL pgbegin(0,plotmode,1,1)
      CALL pgscf(2)
      CALL pgsch(0.7)
      CALL pgask(.FALSE.)
 0100 CONTINUE
      DO i=1,nc
         c(i)=grmin+(grmax-grmin)*REAL(i-1)/REAL(nc-1)
      ENDDO
      CALL bar
      WRITE(*,*)'Object:     ',title
      CALL bar
      WRITE(*,*)'Index - Time(h)'
      column=5
      step=INT(nspec/column)+1
      DO i=1,INT(nspec/column)+1
         WRITE(*,0600)('|',i+j*step,time(i+j*step),j=0,4)
 0600    FORMAT(5(A,I4,F6.1))
      ENDDO
      CALL bar
      WRITE(*,*)'Options: '
      WRITE(*,*)'1 - Continue'
      WRITE(*,*)'2 - Plot Package: ',vpp(pp)
      i=INDEX(wlname,' ')
      WRITE(*,*)'3 - Vel Limits:',vmin,vmax
      WRITE(*,*)'4 - Display Mode: ',plotmode
      WRITE(*,*)'5 - Mapping:      ',vps(ps)
      WRITE(*,*)'E - Exit'
      WRITE(*,*)'>'
      READ(*,5000)opt
      IF (opt.EQ.'1') THEN
         IF (pp.EQ.1) THEN
            CALL main(nspec,n,numdata,first,last,vmin,vmax,time,
     +           title,nline,wl,wlname,pgname,spec,c,nc,ps)
         ELSE
            CALL corrplot(nspec,numdata,vmin,vmax,title,nline,wlname,
     +           pgname,spec,c,nc,ps)
         ENDIF
      ENDIF
      IF (opt.EQ.'2') THEN
         pp=pp+1
         IF (pp.EQ.3) pp=1
      ENDIF
      IF (opt.EQ.'3') THEN
         WRITE(*,*)'Enter vmin,vmax'
         READ(*,*)vmin,vmax
      ENDIF
      IF (opt.EQ.'4') THEN
         CALL pgend
         IF (plotmode.EQ.'/XW')THEN
            plotmode='/PS'
         ELSE
            plotmode='/XW'
         ENDIF
         CALL pgbegin(0,plotmode,1,1)
         CALL pgask(.FALSE.)
      ENDIF
      IF (opt.EQ.'5') THEN
         WRITE(*,*)'Options:'
         DO i=1,3
            WRITE(*,*)i,' - ',vps(i)
         ENDDO
         WRITE(*,*)'Style:'
         READ(*,*)ps
         WRITE(*,*)'Number of levels:'
         READ(*,*)nc
         WRITE(*,*)'Min(0-1):'
         READ(*,*)grmin
         WRITE(*,*)'Max(0-1):'
         READ(*,*)grmax
      ENDIF
      IF (opt.NE.'E') GOTO 0100
 5000 FORMAT(A$)
      END



      SUBROUTINE readparam
      IMPLICIT NONE
      INTEGER nspec
      INTEGER i,nline,nspec
      REAL wl(15),time(150)
      CHARACTER*16 title,wlname(15),pgname(15)
      COMMON /param/ title,nline,nspec,wlname,wl,pgname,time
      CALL find('title:')
      READ(50,5000)title
      CALL find('nline:')
      READ(50,*)nline
      CALL find('nspec:')
      READ(50,*)nspec
      CALL find('wlname:')
      DO i=1,15
         READ(50,5000,ERR=0100)wlname(i)
      ENDDO
 0100 CONTINUE
      CALL find('wl:')
      DO i=1,15
         READ(50,*,ERR=0200)wl(i)
      ENDDO
 0200 CONTINUE
      CALL find('pgname:')
      DO i=1,15
         READ(50,5000,ERR=0300)pgname(i)
      ENDDO
 0300 CONTINUE
 5000 FORMAT(A$)
      RETURN
      END



      SUBROUTINE readline(x,y,n,nspec,numdata,spec)
      IMPLICIT NONE
      CHARACTER*64 spec(15,150)
      REAL x(150,1500),y(150,1500)
      INTEGER numdata,nspec,i,j,n
      DO i=1,nspec
         OPEN(UNIT=50,FILE=spec(n,i))
         numdata=0
         DO j=1,1500
            READ(50,*,END=0100)x(i,j),y(i,j)
            numdata=numdata+1
         ENDDO
 0100    CONTINUE
         CLOSE(50)
      ENDDO
      RETURN
      END



      SUBROUTINE find(word)
      IMPLICIT NONE
      CHARACTER*16 word,fword
      INTEGER i,length
      length=INDEX(word,':')
      i=0
 0100 CONTINUE
      REWIND(50)
 0200 CONTINUE
      i=i+1
      READ(50,5000,END=0100)fword
      IF (fword(1:length).NE.word(1:length))GOTO 0200
 5000 FORMAT(A16)
      RETURN
      END
      

      SUBROUTINE incspec(n,first,last,tvec,time)
      IMPLICIT NONE
      INTEGER first(150),last(150)
      REAL tvec(150),time(150),t0
      INTEGER i,n,start,stop
      CHARACTER opt
      n=0
 0100 CONTINUE
      WRITE(*,*)'First,last,(m = mean, n = individual)(Q to quit):'
      READ(*,*,ERR=0200)start
      READ(*,*,ERR=0200)stop
      READ(*,5000)opt
      IF (opt.EQ.'m') THEN
         n=n+1
         tvec(n)=0.5*(time(start)+time(stop))
         first(n)=start
         last(n)=stop
      ELSE
         DO i=0,(stop-start)
            n=n+1
            tvec(n)=time(start+i)
            first(n)=start+i
            last(n)=start+i
         ENDDO
      ENDIF
      GOTO 0100
 0200 CONTINUE
      t0=tvec(1)
      DO i=1,n
         tvec(i)=tvec(i)-t0
      ENDDO
 5000 FORMAT(A$)
      RETURN
      END



      SUBROUTINE compute(xm,ym,vm,x,y,nspec,n,numdata,nflag,sflag,bias,
     +     ew,snr,first,last,vc,ic)
      IMPLICIT NONE
      REAL ew(150),subspec(1500),bias
      REAL x(150,1500),y(150,1500)
      REAL xm(150,1500),ym(150,1500),vm(150,1500),mm(150,1500)
      REAL imin,imax,vmin,vmax,vc,ic
      INTEGER first(150),last(150)
      INTEGER numdata,nspec,i,j,k,n
      INTEGER nflag,sflag,snr
      DO i=1,nspec
         DO j=1,numdata
            xm(i,j)=x(i,j)+vc
            ym(i,j)=y(i,j)+ic
            IF (nflag.EQ.1)ym(i,j)=ym(i,j)/ew(i)
         ENDDO
      ENDDO
      DO i=1,n
         DO j=1,numdata
            mm(i,j)=0
            DO k=first(i),last(i)
               mm(i,j)=mm(i,j)+ym(k,j)/REAL(last(i)-first(i)+1)
            ENDDO
            vm(i,j)=0
            DO k=first(i),last(i)
               vm(i,j)=vm(i,j)+(ym(k,j)-mm(i,j))**2.
            ENDDO
            vm(i,j)=vm(i,j)/REAL(last(i)-first(i)+1)**2.
            vm(i,j)=vm(i,j)/mm(i,j)
         ENDDO
      ENDDO
      DO j=1,numdata
         subspec(j)=mm(snr,j)
      ENDDO
      k=0
      DO i=1,n
         IF (i.NE.snr)THEN
            k=k+1
            DO j=1,numdata
               ym(k,j)=mm(i,j)-sflag*subspec(j)
               ym(k,j)=ym(i,j)+bias*(k-1)
            ENDDO
         ENDIF
      ENDDO
      n=n-sflag
      CALL max(imax,ym,n,numdata)
      CALL max(vmax,vm,n,numdata)
      CALL min(imin,ym,n,numdata)
      CALL min(vmin,vm,n,numdata)
      IF (imin.GT.imax) imax=ABS(imin)
      IF (vmin.GT.vmax) vmax=ABS(vmin)
      DO i=1,n
         DO j=1,numdata
            ym(i,j)=ym(i,j)/imax
            vm(i,j)=vm(i,j)/vmax
         ENDDO
      ENDDO
      RETURN
      END



      SUBROUTINE main(nspec,n,numdata,first,last,vmin,vmax,time,
     +     title,nline,wl,wlname,pgname,spec,c,nc,ps)
      IMPLICIT NONE
      REAL x(150,1500),y(150,1500),time(150)
      REAL tvec(150),ew(150),c(10),bias,vc,ic,wl(15)
      REAL vmin,vmax,xm(150,1500),ym(150,1500),vm(150,1500)
      INTEGER numdata,first(150),last(150)
      INTEGER sflag,nflag,snr,n,pf,nc,ps,nspec,nline,w,i
      CHARACTER opt
      CHARACTER*3 yn(0:1)
      CHARACTER*64 head,spec(15,150)
      CHARACTER*16 vpf(4),title,wlname(15),pgname(15)
      LOGICAL plex
      w=0
      pf=1
      vc=0
      ic=0
      nflag=0
      sflag=0
      bias=0.
      snr=0
      yn(1)='Yes'
      yn(0)='No'
      vpf(1)='Normal'
      vpf(2)='Variance'
      vpf(3)='Contour'
      vpf(4)='Wtm'
      plex=.FALSE.
 0100 CONTINUE
      CALL bar
      WRITE(*,*)'Options:'
      WRITE(*,*)'1 - Plot'
      WRITE(*,*)'2 - Plot Function:     ',vpf(pf)
      WRITE(*,*)'3 - Line:              ',wlname(w)
      WRITE(*,*)'4 - Included spectra:',n
      WRITE(*,*)'5 - Normalize:         ',yn(nflag)
      WRITE(*,*)'6 - Subtract Spectra:  ',yn(sflag),'(',snr,')'
      WRITE(*,*)'7 - Bias:            ',bias
      WRITE(*,*)'8 - Velocity Shift:  ',vc
      WRITE(*,*)'9 - Intensity Shift: ',ic
      IF (plex) THEN
         WRITE(*,*)'A - Cursor'
      ENDIF
      IF (plex.AND.pf.EQ.1)THEN
         WRITE(*,*)'B - Create ascii files'
      ENDIF
      WRITE(*,*)'M - Up'
      READ(*,5000)opt
      IF (opt.EQ.'2') THEN
         pf=pf+1
         IF (pf.EQ.5)pf=1
      ENDIF
      IF (opt.EQ.'3') THEN
         WRITE(*,*)'Options:'
         DO i=1,nline
            WRITE(*,*)i,' - ',wlname(i)
         ENDDO
         READ(*,*)w
         CALL readline(x,y,w,nspec,numdata,spec)
         head=title(1:INDEX(title,' '))//pgname(w)
      ENDIF
      IF (opt.EQ.'4') CALL incspec(n,first,last,tvec,time)
      IF (opt.EQ.'5') THEN
         nflag=1-nflag
         IF (nflag.EQ.1) CALL equiv(x,y,nspec,numdata,ew,wl)
      ENDIF
      IF (opt.EQ.'6') THEN
         sflag=1-sflag
         snr=0
         IF (sflag.EQ.1) THEN
            WRITE(*,*)'Index of spectra to subtract with:'
            READ(*,*)snr
         ENDIF
      ENDIF
      IF (opt.EQ.'7') THEN
         WRITE(*,*)'Bias: '
         READ(*,*)bias
      ENDIF
      IF (opt.EQ.'8') THEN
         WRITE(*,*)'Velocity Shift(km/s):'
         READ(*,*)vc
      ENDIF
      IF (opt.EQ.'9') THEN
         WRITE(*,*)'Intensity Shift(km/s):'
         READ(*,*)ic
      ENDIF
      IF (opt.EQ.'A') CALL cursor()
      IF (opt.EQ.'B') CALL filew(xm,ym,n,numdata)
      IF (opt.EQ.'M') RETURN
      IF (opt.NE.'1') GOTO 0100
      CALL compute(xm,ym,vm,x,y,nspec,n,numdata,nflag,sflag,bias,
     +     ew,snr,first,last,vc,ic)
      IF (pf.EQ.1) THEN
         CALL plot(xm,ym,n,numdata,vmin,vmax,head)
      ENDIF
      IF (pf.EQ.2) THEN
         CALL variance(xm,vm,n,numdata,vmin,vmax,head)
      ENDIF
      IF (pf.EQ.3) THEN
         CALL interplot(xm,ym,n,numdata,vmin,vmax,tvec,head,c,nc,ps)
      ENDIF
      IF (pf.EQ.4) THEN
         CALL wtm(xm,ym,n,numdata,vmin,vmax,ew,tvec,head,c,nc,ps)
      ENDIF
      plex=.TRUE.
      n=n+sflag
      GOTO 0100
 5000 FORMAT(A$)
      END

      SUBROUTINE plot(xm,ym,n,numdata,vmin,vmax,head)
      REAL xm(150,1500),ym(150,1500)
      REAL x(1500),y(1500),xc,yc
      INTEGER n,numdata
      REAL vmin,vmax,imin,imax
      CHARACTER*64 head
      CALL min(imin,ym,n,numdata)
      CALL max(imax,ym,n,numdata)
      CALL pgsci(5)
      CALL pgenv(vmin,vmax,imin,imax,0,0)
      CALL pglabel('v(km/s)','Intensity',head)
      CALL pgsci(3)
      DO i=1,n
         IF (n.GE.2.AND.n.LE.8) THEN
            CALL pgsci(i)
            xc=vmin+0.02*(vmax-vmin)
            yc=imax-0.25*(imax-imin)*REAL(i)/REAL(n)
            CALL pgpoint(1,xc,yc,ICHAR('o'))
         ENDIF
         DO j=1,numdata
            x(j)=xm(i,j)
            y(j)=ym(i,j)
         ENDDO
         CALL pgline(numdata,x,y)
      ENDDO
      RETURN
      END



      SUBROUTINE variance(xm,vm,n,numdata,vmin,vmax,head)
      IMPLICIT NONE
      REAL xm(150,1500),vm(150,1500)
      REAL imin,imax,vmin,vmax,xc,yc
      INTEGER numdata,i,j,n
      REAL x(1500),y(1500)
      CHARACTER*64 head
      CALL max(imax,vm,n,numdata)
      CALL min(imin,vm,n,numdata)
      CALL pgsci(5)
      CALL pgenv(vmin,vmax,imin,imax,0,0)
      CALL pglabel('v(km/s)','Variance',head)
      CALL pgsci(3)
      DO i=1,n
         IF (n.GE.2.AND.n.LE.8) THEN
            CALL pgsci(i)
            xc=vmin+0.02*(vmax-vmin)
            yc=imax-0.25*(imax-imin)*REAL(i)/REAL(n)
            CALL pgpoint(1,xc,yc,ICHAR('o'))
         ENDIF
         DO j=1,numdata
            x(j)=xm(i,j)
            y(j)=vm(i,j)
         ENDDO
         CALL pgline(numdata,x,y)
      ENDDO
      RETURN
      END

      SUBROUTINE nfind(nmin,nmax,xm,vmin,vmax)
      REAL xm(150,1500),vmin,vmax,vstep
      INTEGER nmin,nmax
      vstep=xm(1,2)-xm(1,1)
      nmin=(vmin-xm(1,1))/vstep+1
      nmax=(vmax-xm(1,1))/vstep+1
      RETURN
      END



      SUBROUTINE wtm(xm,ym,n,numdata,vmin,vmax,ew,tvec,head,c,nc,ps)
      IMPLICIT NONE
      REAL tvec(150)
      REAL xm(150,1500),ym(150,1500),cm(150,1500)
      REAL tr(6),c(10),ew(150),tm(150,150)
      REAL tmax,tmin,maxtime
      REAL vmin,vmax,imax
      INTEGER numdata,i,j,k,n,it,size,nc,ps
      INTEGER nt(100),nmin,nmax
      CHARACTER*64 head
      maxtime=12.
      tmax=0.
      tmin=maxtime
      size=30
      DO i=1,n
         DO j=1,n
            tm(i,j)=0.
         ENDDO
         DO j=1,numdata
            cm(i,j)=0.
         ENDDO
      ENDDO
      DO i=1,n-1
         DO k=i+1,n
            tm(i,k)=tvec(k)-tvec(i)
            IF (tm(i,k).GE.maxtime)tm(i,k)=0.
            IF (tm(i,k).GT.0.) THEN
               IF (tm(i,k).GT.tmax) tmax=tm(i,k)
               IF (tm(i,k).LT.tmin) tmin=tm(i,k)
            ENDIF
         ENDDO
      ENDDO
      imax=0.
      DO j=1,numdata
         DO i=1,100
            nt(i)=0
         ENDDO
         DO i=1,n-1
            DO k=i+1,n
               IF (tm(i,k).GT.0.) THEN
                  it=INT(REAL(size-1)*(tm(i,k)-tmin)/(tmax-tmin)+1.)
                  cm(it,j)=cm(it,j)+ABS(ym(k,j)-ym(i,j))
                  nt(it)=nt(it)+1
               ENDIF
            ENDDO
         ENDDO
         DO i=1,size
            IF (nt(i).GE.1.) THEN
               cm(i,j)=cm(i,j)/nt(i)
            ENDIF
         ENDDO
      ENDDO
      CALL nfind(nmin,nmax,xm,vmin,vmax)
      tr(1)=tmin
      tr(2)=(tmax-tmin)/REAL(size-1)
      tr(3)=0.
      tr(4)=vmin
      tr(5)=0.
      tr(6)=(vmax-vmin)/REAL(numdata-1)
      CALL pgsci(5)
      CALL pgenv(tmin,tmax,vmin,vmax,0,0)
      CALL pglabel('\gDt(h)','v(km/s)',head)
      CALL pgsci(1)
      CALL contourplot(cm,n,numdata,nmin,nmax,c,nc,tr,ps)
 5000 FORMAT(A$)
      RETURN
      END



      SUBROUTINE equiv(x,y,nspec,numdata,ew,wl)
      IMPLICIT NONE
      REAL x(150,1500),y(150,1500)
      REAL vmin,vmax,ew(150),wl,wstep
      INTEGER i,j,numdata,nspec
      wstep=2./3.e5*wl
      WRITE(*,*)'Normalisation vmin, vmax:'
      READ(*,*)vmin,vmax
      DO i=1,nspec
         ew(i)=0.
         DO j=1,numdata
            IF (x(i,j).GE.vmin.AND.x(i,j).LE.vmax) THEN
               ew(i)=ew(i)+(1-y(i,j))*wstep
            ENDIF
         ENDDO
      ENDDO
 5000 FORMAT(A$)
      RETURN
      END



      SUBROUTINE fname(name,i,prefix)
      IMPLICIT NONE
      CHARACTER*4 prefix
      CHARACTER*16 name
      INTEGER i,j,k,l
      j=INT(REAL(i)/10.)
      k=INT(REAL(j)/10.)
      l=INT(REAL(k)/10.)
      name=prefix
     +     //CHAR(MOD(l,10)+48)//char(mod(k,10)+48)
     +     //CHAR(MOD(j,10)+48)//char(mod(i,10)+48)//'.asc'
      RETURN
      END



      SUBROUTINE filew(xm,ym,n,numdata)
      IMPLICIT NONE
      REAL xm(150,1500),ym(150,1500)
      CHARACTER*4 prefix
      CHARACTER*16 filename
      INTEGER i,j,n,numdata
      WRITE(*,*)'Fileprefix(4 letters):'
      READ(*,5000)prefix
      DO i=1,n
         CALL fname(filename,i,prefix)
         OPEN(UNIT=50,FILE=filename)
         DO j=1,numdata
            WRITE(50,*)xm(i,j),ym(i,j)
         ENDDO
         CLOSE(50)
         WRITE(*,*)'Created ',filename
      ENDDO
 5000 FORMAT(A$)
      RETURN
      END



      SUBROUTINE cursor()
      IMPLICIT NONE
      REAL x,y
      CHARACTER ch
      WRITE(*,*)'Press <space> to exit.'
      WRITE(*,*)'     x     y'
 0100 CONTINUE
      CALL pgsci(2)
      CALL pgcurs(x,y,ch)
      CALL pgpoint(1,x,y,ICHAR('0'))
      WRITE(*,0200)x,y
 0200 FORMAT(2(F6.1))
      IF (ch.NE.' ')GOTO 0100
      RETURN
      END



      SUBROUTINE max(tmp,ym,n,numdata)
      IMPLICIT NONE
      REAL ym(150,1500),tmp
      INTEGER numdata,i,j,n
      DO i=1,n
         IF (i.EQ.1)tmp=ym(i,1)
         DO j=1,numdata
            IF (ym(i,j).GT.tmp) tmp=ym(i,j)
         ENDDO
      ENDDO
      RETURN
      END



      SUBROUTINE min(tmp,ym,n,numdata)
      IMPLICIT NONE
      REAL ym(150,1500),tmp
      INTEGER numdata,i,j,n
      DO i=1,n
         IF (i.EQ.1)tmp=ym(i,1)
         DO j=1,numdata
            IF (ym(i,j).LT.tmp) tmp=ym(i,j)
         ENDDO
      ENDDO
      RETURN
      END



      SUBROUTINE contourplot(mx,n,numdata,nmin,nmax,c,nc,tr,ps)
      REAL mx(150,1500),c(10),tr(6),imax
      INTEGER i,j,n,nmin,nmax,nc,ps,numdata
      CALL max(imax,mx,n,numdata)
      DO i=1,n
         DO j=1,numdata
            mx(i,j)=mx(i,j)/imax
         ENDDO
      ENDDO
      IF (ps.EQ.2.OR.ps.EQ.3) THEN
         CALL pggray(mx,150,1500,1,n,nmin,nmax,c(nc),c(1),tr)
      ENDIF
      IF (ps.EQ.1.OR.ps.EQ.3) THEN
         CALL pgcons(mx,150,1500,1,n,nmin,nmax,c,nc,tr)
      ENDIF
      RETURN
      END



      SUBROUTINE interplot(xm,ym,n,numdata,vmin,vmax,tvec,head,c,nc,ps)
      IMPLICIT NONE
      REAL ivec(150),tvec(150),xm(150,1500),ym(150,1500)
      REAL cm(150,1500),tvec(150),tr(6),c(10),vmax,vmin
      INTEGER i,j,n,ns,numdata,iflag,nmax,ps,nc,nmin,nmax
      CHARACTER*64 head
      CHARACTER*3 yn(0:1)
      CHARACTER opt
      yn(0)='No'
      yn(1)='Yes'
      iflag=0
 0100 CONTINUE
      CALL bar
      WRITE(*,*)'Options:'
      WRITE(*,*)'1 - Plot'
      WRITE(*,*)'2 - Interpolate: ',yn(iflag)
      READ(*,0200)opt
 0200 FORMAT(A$)
      IF (opt.EQ.'2') iflag=1-iflag
      IF (opt.EQ.'M') RETURN
      IF (opt.NE.'1') GOTO 0100
      CALL nfind(nmin,nmax,xm,vmin,vmax)
      tr(1)=tvec(1)
      tr(2)=(tvec(n)-tvec(1))/REAL(n-1)
      tr(3)=0.
      tr(4)=vmin
      tr(5)=0.
      tr(6)=(vmax-vmin)/REAL(numdata-1)
      CALL pgsci(5)
      CALL pgenv(tvec(1),tvec(n),vmin,vmax,0,0)
      CALL pglabel('t(h)','v(km/s)',head)
      CALL pgsci(1)
      IF (iflag.EQ.1) THEN
         CALL interp(cm,ym,tvec,ivec,n)
         ns=150
      ELSE
         DO i=1,n
            DO j=1,numdata
               cm(i,j)=ym(i,j)
            ENDDO
         ENDDO
         ns=n
      ENDIF
      CALL contourplot(cm,ns,numdata,nmin,nmax,c,nc,tr,ps)
      RETURN
      END



      SUBROUTINE interp(cm,ym,tvec,ivec,n)
      IMPLICIT NONE
      REAL cm(150,1500),ym(150,1500)
      REAL tvec(n),ivec(n)
      INTEGER i,j,n
      REAL t,iv
      REAL ibiss(9),iprim_0,iprim_n
      iprim_0=0
      iprim_n=0
      DO j=1,1500
         DO i=1,n
            ivec(i)=ym(i,j)
         ENDDO
         CALL spline(tvec,ivec,n,iprim_0,iprim_n,ibiss)
         DO i=1,150
            t=tvec(1)+(tvec(n)-tvec(1))*(i-1)/(150-1)
            CALL splint(tvec,ivec,ibiss,n,t,iv)
            cm(i,j)=iv
         ENDDO
      ENDDO
      RETURN
      END



      SUBROUTINE spline(x,y,n,yp1,ypn,y2)
      INTEGER n,nmax
      REAL yp1,ypn,x(n),y(n),y2(n)
      PARAMETER (nmax=100)
      INTEGER i,k
      REAL p,qn,sig,un,u(nmax)
      IF (yp1.GT..99e30) THEN
         y2(1)=0.
         u(1)=0.
      ELSE
         y2(1)=-0.5
         u(1)=(3./(x(2)-x(1)))*((y(2)-y(1))/(x(2)-x(1))-yp1)
      ENDIF
      DO 0100 i=2,n-1
         sig=(x(i)-x(i-1))/(x(i+1)-x(i-1))
         p=sig*y2(i-1)+2.
         y2(i)=(sig-1.)/P
         u(i)=(6.*((y(i+1)-y(i))/(x(i+
     *        1)-x(i))-(y(i)-y(i-1))/(x(i)-x(i-1)))/(x(i+1)-x(i-1))-sig*
     *        u(i-1))/P
 0100 CONTINUE
      IF (ypn.GT..99e30) THEN
         qn=0.
         un=0.
      ELSE
         qn=0.5
         un=(3./(x(n)-x(n-1)))*(ypn-(y(n)-y(n-1))/(x(n)-x(n-1)))
      ENDIF
      y2(n)=(un-qn*u(n-1))/(qn*y2(n-1)+1.)
      DO 0200 k=n-1,1,-1
         y2(k)=y2(k)*y2(k+1)+u(k)
 0200 CONTINUE
      RETURN
      END



      SUBROUTINE splint(xa,ya,y2a,n,x,y)
      INTEGER N
      REAL x,y,xa(n),y2a(n),ya(n)
      INTEGER k,khi,klo
      REAL a,b,h
      klo=1
      khi=N
 0100 IF (khi-klo.GT.1) THEN
         k=(khi+klo)/2
         IF(xa(k).GT.x)THEN
            khi=k
         ELSE
            klo=k
         ENDIF
         GOTO 0100
      ENDIF
      h=xa(khi)-xa(klo)
      IF (h.EQ.0.) PAUSE 'bad xa input in splint'
      a=(xa(khi)-x)/h
      b=(x-xa(klo))/h
      y=a*ya(klo)+b*ya(khi)+((a**3-a)*y2a(klo)+(b**3-b)*y2a(khi))*(h**
     *     2)/6.
      RETURN
      END



      SUBROUTINE corrplot(nspec,numdata,vmin,vmax,title,nline,wlname,
     +     pgname,spec,c,nc,ps)
      IMPLICIT NONE
      REAL x1(150,1500),y1(150,1500)
      REAL x2(150,1500),y2(150,1500)
      REAL tr(6),c(10)
      REAL xv(150),yv(150),cm(1500,1500)
      REAL r,prob,prob2,z
      REAL vmin,vmax,anti
      INTEGER i,j,k,nspec,numdata,ps,nc,nline
      INTEGER nminx,nmaxx,nminy,nmaxy,xo,yo
      CHARACTER*16 pgname(15),wlname(15),title,vanti(-1:1)
      CHARACTER*64 spec(15,150)
      CHARACTER opt
      LOGICAL plex
      anti=1
      plex=.FALSE.
      vanti(-1)='Negative'
      vanti(+1)='Positive'
      xo=0
      yo=0
 0100 CONTINUE
      CALL bar
      WRITE(*,*)'Options:'
      WRITE(*,*)'1 - Plot'
      WRITE(*,*)'2 - x, y-line:    ',wlname(xo),wlname(yo)
      WRITE(*,*)'3 - Correlation:  ',vanti(anti)
      IF (plex) WRITE(*,*)'4 - Cursor'
      WRITE(*,*)'M - Up'
      READ(*,5000)opt
      IF (opt.EQ.'2') THEN
         DO i=1,nline
            WRITE(*,*)i,' - ',wlname(i)
         ENDDO
         WRITE(*,*)'x-line, y-line:'
         READ(*,*)xo,yo
      ENDIF
      IF (opt.EQ.'3') anti=-1*anti
      IF (opt.EQ.'4')CALL cursor()
      IF (opt.EQ.'M') RETURN
      IF (opt.NE.'1') GOTO 0100
      CALL readline(x1,y1,xo,nspec,numdata,spec)
      CALL readline(x2,y2,yo,nspec,numdata,spec)
      CALL nfind(nminx,nmaxx,x1,vmin,vmax)
      CALL nfind(nminy,nmaxy,x2,vmin,vmax)
      DO i=nminx,nmaxx
         DO k=1,nspec
            xv(k)=y1(k,i)
         ENDDO
         DO j=nminy,nmaxy
            DO k=1,nspec
               yv(k)=y2(k,j)
            ENDDO
            r=0.
            IF (xv(1).NE.0..OR.yv(1).NE.0.) THEN
               CALL pearsn(xv,yv,nspec,r,prob,prob2,z)
            ENDIF
            cm(i,j)=anti*r
         ENDDO
      ENDDO
      tr(1)=vmin
      tr(2)=(vmax-vmin)/REAL(numdata-1)
      tr(3)=0.
      tr(4)=vmin
      tr(5)=0.
      tr(6)=(vmax-vmin)/REAL(numdata-1)
      CALL pgsci(5)
      CALL pgenv(vmin,vmax,vmin,vmax,1,-2)
      CALL pgbox('BCNTS',100.,2,'BCNTS',100.,2)
      CALL pglabel(pgname(xo),pgname(yo),title)
      CALL pgsci(1)
      IF (ps.EQ.2.OR.ps.EQ.3) THEN
         CALL pggray(cm,1500,1500,nminx,nmaxx,nminy,nmaxy,c(nc),c(1),tr)
      ENDIF
      IF (ps.EQ.1.OR.ps.EQ.3) THEN
         CALL pgcons(cm,1500,1500,nminx,nmaxx,nminy,nmaxy,c,nc,tr)
      ENDIF
      plex=.TRUE.
      GOTO 0100
 5000 FORMAT(A$)
      END



      SUBROUTINE pearsn(x,y,n,r,prob,prob2,z)
      PARAMETER (tiny=1.e-20)
      DIMENSION x(n),y(n)
      ax=0.
      ay=0.
      DO 0100 j=1,n
         ax=ax+x(j)
         ay=ay+y(j)
 0100 CONTINUE
      ax=ax/n
      ay=ay/n
      sxx=0.
      syy=0.
      sxy=0.
      DO 0200 j=1,n
         xt=x(j)-ax
         yt=y(j)-ay
         sxx=sxx+xt**2
         syy=syy+yt**2
         sxy=sxy+xt*yt
 0200 CONTINUE
      r=sxy/SQRT(sxx*syy)
      z=0.5*alog(((1.+r)+tiny)/((1.-r)+tiny))
      df=n-2
      t=r*SQRT(df/(((1.-r)+tiny)*((1.+r)+tiny)))
      prob=0.1                  !betai(0.5*df,0.5,df/(df+t**2))
      prob2=0.1                 !erfcc(ABS(z*SQRT(n-1.))/1.4142136)
      RETURN
      END



      FUNCTION betai(a,b,x)
      IF(x.LT.0..or.x.GT.1.)PAUSE 'bad argument X in BETAI'
      IF(x.EQ.0..or.x.eq.1.)THEN
         bt=0.
      ELSE
         bt=EXP(gammln(a+b)-gammln(a)-gammln(b)
     *        +a*alog(x)+b*alog(1.-x))
      ENDIF
      IF(x.LT.(a+1.)/(a+b+2.))THEN
         betai=bt*betacf(a,b,x)/A
         RETURN
      ELSE
         betai=1.-bt*betacf(b,a,1.-x)/B
         RETURN
      ENDIF
      END




      FUNCTION erfcc(x)
      z=ABS(x)
      t=1./(1.+0.5*z)
      erfcc=t*EXP(-z*z-1.26551223+t*(1.00002368+t*(.37409196+
     *     t*(.09678418+t*(-.18628806+t*(.27886807+t*(-1.13520398+
     *     t*(1.48851587+t*(-.82215223+t*.17087277)))))))))
      IF (x.LT.0.) erfcc=2.-erfcc
      RETURN
      END



      FUNCTION betacf(a,b,x)
      PARAMETER (itmax=100,eps=3.e-7)
      am=1.
      bm=1.
      az=1.
      qab=a+B
      qap=a+1.
      qam=a-1.
      bz=1.-qab*x/qap
      DO 0100 m=1,itmax
         em=M
         tem=em+em
         d=em*(b-m)*x/((qam+tem)*(a+tem))
         ap=az+d*am
         bp=bz+d*bm
         d=-(a+em)*(qab+em)*x/((a+tem)*(qap+tem))
         app=ap+d*az
         bpp=bp+d*bz
         aold=az
         am=ap/bpp
         bm=bp/bpp
         az=app/bpp
         bz=1.
         IF(ABS(az-aold).LT.eps*abs(az))GO TO 0200
 0100 CONTINUE
      PAUSE 'A or B too big, or ITMAX too small'
 0200 betacf=az
      RETURN
      END



      FUNCTION gammln(xx)
      REAL*8 cof(6),stp,half,one,fpf,x,tmp,ser
      DATA cof,stp/76.18009173d0,-86.50532033d0,24.01409822d0,
     *     -1.231739516d0,.120858003d-2,-.536382d-5,2.50662827465d0/
      DATA half,one,fpf/0.5d0,1.0d0,5.5d0/
      x=xx-one
      tmp=x+fpf
      tmp=(x+half)*LOG(tmp)-tmp
      ser=one
      DO 0100 j=1,6
         x=x+one
         ser=ser+cof(j)/x
 0100 CONTINUE
      gammln=tmp+LOG(stp*ser)
      RETURN
      END

      SUBROUTINE bar
      WRITE(*,*)'=============================='
      RETURN
      END
