Fortran 90 Toolkit to Accompany "A Class of Fast Methods for
      Processing Irregularly Sampled or Otherwise Inhomogeneous
      One-Dimensional Data" by George B. Rybicki and William H. Press

What this file is:

This file contains several Fortran 90 modules, with *very* brief
explanations, which collectively implement the methods described in
the above-named paper.  The purpose of this file is to make
algorithmically precise the necessarily abbreviated descriptions in
the paper.  It is intended that the user treat this file as a
"toolkit", using bits and pieces as needed for implementing his/her
own applications.

What this file IS NOT:

It is NOT polished, ready-to-use software.  It is NOT guaranteed to be
free of all bugs, or to treat properly all exceptional cases that
might arise.  It is NOT the authors' last word on this subject -- our
development of algorithms and applications continues.  (This file, as
accompanying the above-named paper, will NOT be updated after the
paper is published in final form.)  And, it is NOT a programming
text in Fortran 90.

Why Fortran 90?

A couple of reasons.  (1) That language's vector and matrix
constructions allow for quite concise (and potentially quite elegant)
expression of our algorithms.  (2) This research happened to overlap
in time one author's (WHP) work a forthcoming new version of the book
"Numerical Recipes", which will be in Fortran 90.  The best
introduction to Fortran 90 is the book "Fortran 90 Explained" by
Metcalf and Reid (Oxford University Press, 1992).  Fortran 90
compilers for most types of workstations are in late stages of
development and will be available (if not already) by late 1994.  This
work was done on a DEC 3000/600 (Alpha) workstation, using a Field
Test version of DEC Fortran 90.

Copyright Notice:

Portions of this file (those derived from the Numerical Recipes books)
are Copyright 1986-1992 by Numerical Recipes Software.  Permission is
granted for the unlimited noncommercial use of the routines included
herein, PROVIDED THAT such use (unless authorized by an additional
license) is exclusively in support of, and incidental to, the use of
the algorithms described in the above-named paper.

********file reestimate.f90*****************************

! This module contains the user-callable routines remap and reest.

! Remap "rectifies" irregularly spaced data (y) known (with errors sig) on
! one set of x's to data (yout) on a different, generally interleaved,
! set of xout's.  The data is characterized by a single scalar population
! variance (psig**2) and a single inverse decorrelation distance w.
! All quantities are input, except yout, which is output.

! Reest is a lower-level routine that reestimates a single data set
! (x,y,sig) onto a set of improved y values yout, located at the same
! data points x.  Now psig and w are general vectors.  The optional
! argument ybar, if present, causes an unbiased (Gauss-Markov) estimate
! to be used (essentially removing and then restoring the mean of the
! data) and returns the value of the data's mean.  The optional argument
! chisq, if present, causes the chi-square of the data's internal
! consistency to be returned as chisq.

! A trivial test routine is appended (commented out).

MODULE reestimate
CONTAINS
	SUBROUTINE remap(x,y,sig,psig,w,xout,yout)
	IMPLICIT NONE
	INTEGER :: n,nout,i,j,k,m,ntot
	REAL :: psig,w,ave
	REAL, DIMENSION(:) :: x,y,sig,xout,yout
	REAL, DIMENSION(:), ALLOCATABLE :: xx,yy,ppsig,ssig,ww,yyout
	INTEGER, DIMENSION(:), ALLOCATABLE :: iflag
	n=size(x)
	nout=size(xout)
	if (size(y) /= n .or. size(sig) /= n .or. size(yout) /= nout) &
		 call nrerror('remap: wrong sizes')
	if (any(x(2:n)<=x(1:n-1)).or.any(xout(2:nout)<=xout(1:nout-1))) &
		call nrerror('remap: values must be distinct and ordered')
	ave=sum(y)/n
	ntot=n+nout
	allocate(xx(ntot),yy(ntot),ppsig(ntot),ssig(ntot),ww(ntot), &
		yyout(ntot),iflag(ntot))
	iflag=0
	j=1
	k=1
	do m=1,ntot
		if(j>n.or.k>nout) exit
		if (x(j)<=xout(k)) then
			iflag(m)=j
			j=j+1
		else ! (x(j)>xout(k))
			iflag(m)=-k
			k=k+1
		endif
	enddo
	if (j<=n) iflag(m:ntot)=(/ (i,i=j,n) /)
	if (k<=nout) iflag(m:ntot)=(/ (-i,i=k,nout) /)
	where (iflag>0)
		xx=x(iflag)
		yy=y(iflag)
		ssig=sig(iflag)
		ppsig=psig
		ww=w
	elsewhere
		xx=xout(-iflag)
		yy=ave
		ssig=100.*psig
		ppsig=psig
		ww=w
	end where
	call reest(xx,yy,ppsig,ssig,ww(1:ntot-1),yyout)
	yout=pack(yyout,iflag<0)
	END SUBROUTINE

	SUBROUTINE reest(x,y,psig,sig,w,yout,ybar,chisq)
	IMPLICIT NONE
	REAL, DIMENSION(:) :: x,y,psig,sig,w,yout
	REAL, OPTIONAL :: ybar,chisq
	REAL, PARAMETER :: TINY=1.e-4
	INTEGER :: n
	REAL, DIMENSION(size(x)) :: b,temp,bb,e
	REAL, DIMENSION(size(x)-1) :: a,c,r,aa,cc
	n=size(x)
	if (size(y) /= n.or.size(psig) /= n.or.size(sig) /= n.or.size(yout) &
		/= n.or.size(w) /= n-1) call nrerror('reest: wrong sizes')
	if (any(x(2:n)<x(1:n-1))) &
		call nrerror('reest: x values must be distinct and ordered')
	temp=sig**2
! compute Rybicki's tridiagonal elements
	r = exp(-abs(w*(x(2:n)-x(1:n-1))))
	a = -1./max(1./r-r,TINY)
	r = r*a
	b(1) = 1.-r(1)
	b(2:n-1) = 1.-r(2:n-1)-r(1:n-2)
	b(n) = 1.-r(n-1)
	a = a/(psig(1:n-1)*psig(2:n))
	b = b/psig**2
! mult by diagonal matrix sig**2 and add identity
	aa = a*temp(2:n)
	bb = b*temp+1.
	cc = a*temp(1:n-1)
	if (present(ybar)) then
! calculate (S+N)^(-1)*E using temp for scratch
		call tridag(aa,bb,cc,spread(1.,1,n),temp)
		e=symtrimul(a,b,temp)
		ybar=dot_product(e,y)/sum(e)
		temp=y-ybar
	else
		temp=y
	endif
! calculate S*(S+N)^(-1)*(y-ybar)
!	do j=1,10
!		write (*,*) j,aa(j),bb(j),cc(j),temp(j)
!	enddo
	call tridag(aa,bb,cc,temp,yout)
	if (present(chisq)) then
		chisq=dot_product(temp,symtrimul(a,b,yout))
	endif
	if (present(ybar)) then
		yout=yout+ybar
	endif
	END SUBROUTINE reest

	FUNCTION symtrimul(a,b,u)
	REAL, DIMENSION(:) :: a,b,u
	REAL, DIMENSION(size(u)) :: symtrimul
	INTEGER :: n
	n=size(u)
	symtrimul(1)=b(1)*u(1)+a(1)*u(2)
	symtrimul(2:n-1)=a(1:n-2)*u(1:n-2)+b(2:n-1)*u(2:n-1)+a(2:n-1)*u(3:n)
	symtrimul(n)=a(n-1)*u(n-1)+b(n)*u(n)
	END FUNCTION symtrimul

	RECURSIVE SUBROUTINE tridag(a,b,c,r,u)
	IMPLICIT NONE
	REAL, DIMENSION(:) :: a,b,c,r,u
	INTEGER, PARAMETER :: NPAR=1000000000
	INTEGER :: n,n2,nm,nx,j
	REAL, DIMENSION(:), ALLOCATABLE :: gam,x,y,z,q,piva,pivc
	REAL :: bet
	n=size(b)
	if (size(a) /= n-1 .or. size(c) /= n-1 .or. size(r) /= n &
		.or. size(u) /= n) call nrerror('bad sizes in tridag')
	if (n<NPAR) then
		bet=b(1)
		if(bet==0.) call nrerror('Error 1 in tridag')
		allocate( gam(n) )
		u(1)=r(1)/bet
		do j=2,n
			gam(j)=c(j-1)/bet
			bet=b(j)-a(j-1)*gam(j)
			if(bet==0.) call nrerror('Error 2 in tridag')
			u(j)=(r(j)-a(j-1)*u(j-1))/bet
		enddo
		do j=n-1,1,-1
			u(j)=u(j)-gam(j+1)*u(j+1)
		enddo
		deallocate(gam)
	else
		call nrerror('this tridag has its parallel leg amputated')
	endif
	return
	END SUBROUTINE tridag

	SUBROUTINE nrerror(string)
	CHARACTER(LEN=*) :: string
	write (*,*) 'nrerror: ',string
	STOP 'program terminated by nrerror'
	END SUBROUTINE nrerror

END MODULE

!	program test_remap
!	USE reestimate
!	REAL, DIMENSION(100) :: x,y,sig
!	REAL, DIMENSION(99) :: xout,yout
!	x=(/ (1.0*i,i=1,100) /)
!	xout=x(1:99)+0.0
!	y=sin(0.5*x)
!	sig=.1
!	psig=1.
!	w=1./12.
!	call remap(x,y,sig,psig,w,xout,yout)
!	do j=1,100
!		write (*,*) x(j),y(j),' in'
!	enddo
!	do j=1,99
!		write (*,*) xout(j),yout(j),' out'
!	enddo
!	END

*************end file******************************************

*************file lopast.f90***********************************

! The function lopast(y,x,f) returns a vector that is the low-pass
! filtered version of the data points y, at the abcissas x, with a 3dB
! cutoff frequency of f.  This is pretty much exactly as described in
! the paper.

! High pass filtering is obtained by changing the line
!	lopast=y-real(ans)
! to
!	lopast=real(ans)
! and (if you want f to be the 3dB point again) changing the
! constant 5.5380744 to 3.56427.  (If the filters were characterized
! by 6dB cutoffs instead of 3dB cutoffs, the constants would of course
! be the same -- this is only semantics in the definition of cutoff.)

! A trivial test routine is appended (commented out).

MODULE lopast_module
CONTAINS
	FUNCTION lopast(y,x,f)
	IMPLICIT NONE
	REAL, DIMENSION(:) :: y,x
	REAL :: f
	REAL, DIMENSION(size(y)) :: lopast
	INTEGER :: n
	COMPLEX :: fac
	COMPLEX, DIMENSION(size(y)) :: b,rhs,ans
	COMPLEX, DIMENSION(size(y)-1) :: r,e,w,re
	n=size(y)
	if (size(x) /= n) pause 'bad sizes in lopast'
	fac=abs(f)*5.5380744*cmplx(1.,1.)
	w=fac*abs(x(2:n)-x(1:n-1))
	r=exp(-w)  ! note this is a complex exponential !
	e=1./(1./r - r)
	re=r*e
	b=1.
	b(1:n-1)=b(1:n-1)+re
	b(2:n)=b(2:n)+re
	w = 0.5/w
	rhs(1:n-1)=w*(y(1:n-1)-y(2:n))
	rhs(n)=0.
	rhs(2:n)=rhs(2:n)+w*(y(2:n)-y(1:n-1))
	e=-e
	call ctridag(e,b,e,rhs,ans)
	lopast=y-real(ans)
	END FUNCTION lopast

	SUBROUTINE ctridag(a,b,c,r,u)
	IMPLICIT NONE
	COMPLEX, DIMENSION(:) :: a,b,c,r,u
	INTEGER :: n,j
	COMPLEX, DIMENSION(size(b)) :: gam
	COMPLEX :: bet
	n=size(b)
	if (size(a) /= n-1 .or. size(c) /= n-1 .or. size(r) /= n &
		.or. size(u) /= n) pause 'bad sizes in ctridag'
	bet=b(1)
	if(bet==(0.,0.)) pause 'Error 1 in ctridag'
	u(1)=r(1)/bet
	do j=2,n
		gam(j)=c(j-1)/bet
		bet=b(j)-a(j-1)*gam(j)
		if(bet==(0.,0.)) pause 'Error 2 in ctridag'
		u(j)=(r(j)-a(j-1)*u(j-1))/bet
	enddo
	do j=n-1,1,-1
		u(j)=u(j)-gam(j+1)*u(j+1)
	enddo
	END SUBROUTINE ctridag
END MODULE

!	PROGRAM test_lopast
!	USE lopast_module
!	INTEGER, PARAMETER :: n=200
!	REAL, DIMENSION(n) :: x,y,a,b,c,d
!	x = (/ (j+0.5*sin(0.1*real(j)),j=1,n) /)
!	y = 0.
!	y(1:20)=1.
!	a=lopast(y,x,0.01)
!	b=lopast(y,x,0.03)
!	c=lopast(y,x,0.1)
!	d=lopast(y,x,0.3)
!	do j=1,n
!		write (*,*) x(j),a(j),b(j),c(j),d(j)
!	enddo
!	END

*******************end file************************

***************file fast_fit.f90*******************

! The routine fastfit returns a vector of parameters q and a chi square
! value chisq that describe the fit of a linear combination of the
! columns of the matrix el (with abscissas xel) to the data set x,y,sig
! (having a different set of abscissas x).  psig and pw are the same as
! psig and w in the previous modules.

! A trivial test routine is appended (commented out).

MODULE fast_fitting
CONTAINS
	SUBROUTINE fastfit(x,y,sig,psig,pw,xel,el,q,chisq)
! note psig and pw are scalars in this version
	REAL :: psig,pw,chisq
	REAL, DIMENSION(:) :: x,y,sig,xel,q
	REAL, DIMENSION(:,:) :: el
	REAL, DIMENSION(size(q),1) :: qq
	REAL, DIMENSION(size(x)+size(xel),size(el,2)+1) :: ell,cil
	REAL, DIMENSION(size(x)+size(xel)) :: xx
	LOGICAL, DIMENSION(size(x)+size(xel)) :: mx
	INTEGER :: n,m,nl,j
	n=size(x)
	nl=size(xel)
	m=size(el,2)
	if (size(y) /= n .or. size(sig) /= n .or. size(el,1) /= nl .or. &
		size(q) /= m) call nrerror('fastfit: bad sizes')
	call mux(x,xel,mx)
	xx=unpack(x,mx,0.)+unpack(xel,.not.mx,0.)
	ell(:,m+1)=unpack(y,mx,0.)
	do j=1,m
		ell(:,j)=unpack(el(:,j),.not.mx,0.)
	enddo
	call apply_cinv(xx,ell,unpack(sig,mx,0.),psig,pw,cil)
	qq(:,1)=matmul(transpose(ell(:,1:m)),cil(:,m+1))
	call gaussj( matmul(transpose(ell(:,1:m)),cil(:,1:m)) , qq)
	q=-qq(:,1)
	ell(:,1)=ell(:,m+1)+matmul(ell(:,1:m),q)
	cil(:,1)=cil(:,m+1)+matmul(cil(:,1:m),q)
	chisq=dot_product(ell(:,1),cil(:,1))
	END SUBROUTINE

	SUBROUTINE apply_cinv(x,y,sig,psig,pw,yout)
! note psig and pw are scalars in this version
	IMPLICIT NONE
	REAL :: psig,pw
	REAL, DIMENSION(:) :: x,sig
	REAL, DIMENSION(:,:) :: y,yout
	REAL, PARAMETER :: TINY=1.e-4
	INTEGER :: n,m,j
	REAL, DIMENSION(size(x)) :: b,temp,bb
	REAL, DIMENSION(size(x)-1) :: a,r,aa,cc
	n=size(x)
	m=size(y,2)
	if (size(y,1) /= n .or. size(sig) /= n .or. size(yout,1) /= n .or. &
		size(yout,2) /= m) call nrerror('apcinv: wrong sizes')
	if (any(x(2:n)<x(1:n-1))) &
		call nrerror('apcinv: x values must be distinct and ordered')
	temp=sig**2
! compute elements of the (tridiagonal) inverse
	r = exp(-abs(pw*(x(2:n)-x(1:n-1))))
	a = -1./max(1./r-r,TINY)
	r = r*a
	b(1) = 1.-r(1)
	b(2:n-1) = 1.-r(2:n-1)-r(1:n-2)
	b(n) = 1.-r(n-1)
	a = a/psig**2
	b = b/psig**2
! construct (1+NS^[-1]) into aa,bb,cc
	aa = a*temp(2:n)
	bb = b*temp+1.
	cc = a*temp(1:n-1)
	do j=1,m
! mult y by (1+NS^[-1])^[-1]
		call tridag(aa,bb,cc,y(:,j),temp)
! mult prev result by S^[-1]
		yout(:,j)=symtrimul(a,b,temp)
	enddo
	END SUBROUTINE apply_cinv

	FUNCTION symtrimul(a,b,u)
	REAL, DIMENSION(:) :: a,b,u
	REAL, DIMENSION(size(u)) :: symtrimul
	INTEGER :: n
	n=size(u)
	symtrimul(1)=b(1)*u(1)+a(1)*u(2)
	symtrimul(2:n-1)=a(1:n-2)*u(1:n-2)+b(2:n-1)*u(2:n-1)+a(2:n-1)*u(3:n)
	symtrimul(n)=a(n-1)*u(n-1)+b(n)*u(n)
	END FUNCTION symtrimul

	RECURSIVE SUBROUTINE tridag(a,b,c,r,u)
	IMPLICIT NONE
	REAL, DIMENSION(:) :: a,b,c,r,u
	INTEGER, PARAMETER :: NPAR=1000000000
	INTEGER :: n,j
	REAL, DIMENSION(:), ALLOCATABLE :: gam
	REAL :: bet
	n=size(b)
	if (size(a) /= n-1 .or. size(c) /= n-1 .or. size(r) /= n &
		.or. size(u) /= n) call nrerror('bad sizes in tridag')
	if (n<NPAR) then
		bet=b(1)
		if(bet==0.) call nrerror('Error 1 in tridag')
		allocate( gam(n) )
		u(1)=r(1)/bet
		do j=2,n
			gam(j)=c(j-1)/bet
			bet=b(j)-a(j-1)*gam(j)
			if(bet==0.) call nrerror('Error 2 in tridag')
			u(j)=(r(j)-a(j-1)*u(j-1))/bet
		enddo
		do j=n-1,1,-1
			u(j)=u(j)-gam(j+1)*u(j+1)
		enddo
		deallocate(gam)
	else
		call nrerror('this tridag has its parallel leg amputated')
	endif
	return
	END SUBROUTINE tridag

	SUBROUTINE nrerror(string)
	CHARACTER(LEN=*) :: string
	write (*,*) 'nrerror: ',string
	STOP 'program terminated by nrerror'
	END SUBROUTINE nrerror

	SUBROUTINE gaussj(a,b)
	REAL, DIMENSION(:,:) :: a,b
	LOGICAL, DIMENSION(size(a,1),size(a,1)) :: mask
	LOGICAL, DIMENSION(size(b,1),size(b,2)) :: maskb
	INTEGER, DIMENSION(size(a,1)) :: ipiv,indxr,indxc
	REAL, DIMENSION(size(a,1)) :: dumc
	INTEGER, TARGET :: irc(2)
	INTEGER :: m,n
	INTEGER, POINTER :: irow,icol
	irow => irc(1)
	icol => irc(2)
	n=size(a,1)
	m=size(b,2)
	ipiv=0
	do i=1,n
		mask = spread(ipiv==0,1,n).and.spread(ipiv==0,2,n)
		irc=maxloc(abs(a),mask)
		ipiv(icol)=ipiv(icol)+1
		if(count(ipiv>1)/=0) call nrerror('singular matrix')
		if(irow/=icol) then
			call SWAP(a(irow,:),a(icol,:))
			call SWAP(b(irow,:),b(icol,:))
		endif
		indxr(i)=irow
		indxc(i)=icol
		if(a(icol,icol)==0.) call nrerror('SINGULAR MATRIX')
		pivinv=1./a(icol,icol)
		a(icol,icol)=1.
		a(icol,:)=a(icol,:)*pivinv
		b(icol,:)=b(icol,:)*pivinv
		dumc=a(:,icol)
		a(:,icol)=0.
		a(icol,icol)=pivinv
		mask=.true.
		mask(icol,:)=.false.
		where(mask) a=a-spread(a(icol,:),1,n)*spread(dumc,2,n)
		maskb=.true.
		maskb(icol,:)=.false.
		where(maskb) b=b-spread(b(icol,:),1,n)*spread(dumc,2,m)
	enddo	
	do l=n,1,-1
		call SWAP(a(:,indxr(l)),a(:,indxc(l)))
	enddo
	return
	END SUBROUTINE gaussj

	SUBROUTINE swap(a,b)
	REAL, DIMENSION(:) :: a,b
	REAL, DIMENSION(SIZE(a)) :: dum
	dum=a
	a=b
	b=dum
	END SUBROUTINE swap

	SUBROUTINE mux(x1,x2,flag)
	IMPLICIT NONE
	REAL, DIMENSION(:) :: x1,x2
	LOGICAL, DIMENSION(:) :: flag
	INTEGER :: n1,n2,nt,j,k,m
	n1=size(x1)
	n2=size(x2)
	nt=n1+n2
	if (size(flag) /= nt) call nrerror('mux: flag wrong size')
	if (any(x1(2:n1)<x1(1:n1-1)).or.any(x2(2:n2)<x2(1:n2-1))) &
		call nrerror('mux: x1 and x2 must each be ordered')
	j=1
	k=1
	do m=1,nt
		if(j>n1.or.k>n2) exit
		if (x1(j)<=x2(k)) then
			flag(m)=.true.
			j=j+1
		else
			flag(m)=.false.
			k=k+1
		endif
	enddo
	if (j<=n1) then
		flag(m:nt)=.true.
	else if (k<=n2) then
		flag(m:nt)=.false.
	endif
	END SUBROUTINE mux
END MODULE

!	program test_fast_fitting
!	USE fast_fitting
!	INTEGER, PARAMETER :: N=200,NEL=199
!	REAL :: x(N),y(N),sig(N),xel(NEL),el(NEL,3),q(3)
!1	write (*,*) 'guess f (20 is right):'
!	read (*,*) f
!	eps=0.001
!	x=(/ (i+0.5*sin(197.*i),i=1,N) /)
!	xel=(/ (i+.33,i=1,NEL) /)
!	y=7.0*sin(x/20.)+5.0*sin(x/30.) +eps*sin(1000.*x)
!	el(:,1)=sin(xel/f)
!	el(:,2)=sin(xel/30.)
!	el(:,3)=1.
!	sig=eps
!	psig=100.
!	pw=1./30.
!	call fastfit(x,y,sig,psig,pw,xel,el,q,chisq)
!	write (*,*) 'q=',q,'  chisq=',chisq
!	goto 1
!	END

******************end file***************************