!                    BSPLINES 

! Purpose: evaluate B-spline basis functions
!          linear, quadratic or cubic
!          equidistant (default) or at user-specified knots

! Input: via command line options (-h will list them) or interactively
!        program will prompt for esential options not given on command line;
!        ages and knot positions can be given in separate files

! Choose
!   degree: command line option -d1, -d2 or -d3 for linear, quadratic or cubic;
!           if not given, program will prompt for it            
!   "ages": command line option -a to read min. & max. age and step size 
!           interactively; if not given program will prompt for it,
!           if given, program will expect to read ages from file "agesfile",
!           default name "ages.dat"
!   knots:  command line option -eqx specificies a spline function with equal
!           intervals, defined by x knots (x must be integer); if not given,
!           program will prompt for the no. of knots.
!           to select a spline with user-selected knots, a "knotsfile" (default
!           name "knots" must exist in the working directory; program will
!           count no. of knots in the file

! Format of "agesfile":
!          one age per line
!          each line comprised of two space-separated values: running no. 
!          (integer) and value for age (integer)

! Format for "knotsfile":
!          one line per knot containing a single (real) value,
!          knots should be in ascending order;
!          it is the user's responsibility to ensure that knots are sensible &
!          that the external knots & the range of ages specified are consistent!

! Output file "bspfile" with default name "bspline_Xn.dat" where X is the 
!          degree of the spline (L, Q or C) and n is the no. of knots given:
!          lines in order of ages 
!          one line per age with n coefficients (real), followed by value for
!          age; adjust format for output as necessary

! Compilation note: 'gfortran' may not accommadate the intrinsic functions
!          'iargc' and 'getarg' in subroutine cmdline; replace by
!          'k = command_argument_count()' and 
!          'call get_command_argument( iarg, arg )' respectively

! ======================================================================
  program Bsplines
! ===================================================@C-km-2004-========

  implicit none
! -------------------- change as desired ----------------------------------
! default name for file with user-chosen knots
  character(len=20)                      :: knotsfile ='knots              '
! default name for file with "ages" at which to evaluate B-spline function
  character(len=20)                      :: agesfile = 'ages.dat           '
! default name for output file
  character(len=16)                      :: bsplfile = 'bspline_X00.dat'
! format for output file
  character(len=12)                      :: fmt = '(00f12.8,i8)'
! ------------------- end of changes --------------------------------------
  character(len=1), dimension(3)         :: degree = (/ 'L', 'Q', 'C' /)
  real(8), dimension(:), allocatable     :: kvec, tvec, w, p
  integer, dimension(:), allocatable     :: iiage
  real(8), dimension(:,:,:), allocatable :: mm
  integer                                :: nknot, ndeg, nord, mknot, nage, nrr
  integer                                :: ii, i, j, iage, ia, iage1, iage2, &
&                                           istep
  real(8)                                :: xx, age, age1, age2, aa
  logical                                :: verbose, degset, kntset, equal,  &
&                                           ages, lex

! read command line options
  call cmdline

! degree of basis functions
 if(.not.degset)then
    write(*,*)'degree ? (1=linear, 2=quad., 3=cubic)'
    read(*,*) ndeg; if(ndeg <1 .or. ndeg >3 ) stop 'wrong degree'
 end if
 nord = ndeg+1  ! order of fit

! ages for which to evaluate basis functions
100 continue
if( ages ) then
   write(*,*)'Minimum age ?';   read(*,*)iage1
   write(*,*)'Maximum age ?';   read(*,*)iage2
   if(iage2.le.iage1)then
      write(*,*)'error : max. age must be greater than min. age'
      stop 'Bsplines'
   end if
   write(*,*)'Step size ?';  read(*,*) istep
   allocate( iiage(iage2-iage1+1) ); nage=0
   do iage=iage1,iage2,istep
      nage=nage+1; iiage(nage) = iage
   end do
   if(iiage(nage)<iage2)write(*,*) &
&                       'Warning : Age range not a multiple of step size !'
else
   inquire(file = agesfile,exist=lex)
   if( .not. lex)then
      write(*,*)'File',agesfile,'does not exist !'
      ages=.true.;  go to 100
   end if
   open(1,file=agesfile,status='old', action='read' )
   nage = 0
   do
      read(1,*,iostat=ii) xx;  if(ii /= 0) exit
      nage = nage + 1
   end do
   allocate( iiage(nage));  rewind(1)
   do
      read(1,*,iostat=ii) nage, iiage(nage);   if(ii/=0)exit
      if(verbose) write(*,*)'age', nage, iiage(nage)
   end do
   close(1)
end if
if( verbose ) then
   write(*,'(a,t30,1h=,i8)')'No. of ages',nage
   write(*,'(a,t30,1h=,2i8)')'Age range',iiage(1),iiage(nage)
end if

! where are the knots ?
200 continue
if( equal ) then
   if( .not. kntset ) then
       write(*,'(a)')'No. of knots ?'; read(*,*) nknot; kntset=.true.
   end if
   age1=iiage(1);  age2=iiage(nage); aa= ( age2 - age1 ) / dble( nknot - 1 )
   allocate( kvec(nknot) )
   do i = 1, nknot
      kvec(i) = age1 + (i-1)*aa
   end do
   mknot = nknot + 2*ndeg;  nrr = nknot + ndeg -1
   allocate( tvec(mknot) )
   tvec(nord:nord+nknot-1) = kvec
   do i=1,ndeg
      tvec(i)= age1 - (ndeg-i+1)*aa
   end do
   do i=1,ndeg
      tvec(nknot+nord-1+i) = age2 + i* aa
   end do
else
   inquire(file=knotsfile,exist=lex)
   if(.not.lex) then
      write(*,'(a,1h",a,1h",a)')'File ',trim(knotsfile),' does not exist !'
      equal=.true.; go to 200
   end if
   open(1,file=knotsfile,status='old',action='read') ;  nknot=0
   do
      read(1,*,iostat=ii)xx; if(ii /= 0)exit
      nknot = nknot+1
   end do
   write(*,*) 'no. of user-specified knots =', nknot
   allocate( kvec(nknot) ); rewind(1)
   do i=1,nknot
      read(1,*,iostat=ii)kvec(i);  if(ii /= 0)exit
      if(verbose)write(*,*)'knot',i,kvec(i)
   end do
   close(1)
   if(iiage(1)<kvec(1) .or. iiage(nage) .ge. kvec(nknot) ) then
       ! check that range of ages & external knots are consistent
       write(*,*)'range of ages and knots specified disagree !'
       write(*,*)'external knots given =',kvec(1),kvec(nknot)
       write(*,*)'range of ages found  =',iiage(1),iiage(nage)
       stop 'Bsplines'
   end if
   mknot = nknot + 2*ndeg;  nrr = nknot + ndeg -1  ! no. of RR coefficients
   allocate( tvec(mknot) )
   tvec(1:ndeg)=kvec(1); tvec(nord:nord+nknot-1) = kvec
   tvec(nord+nknot:mknot) = kvec(nknot)
end if
if(verbose)write(*,*)'(Extended) Vector of knots ',tvec

! basis functions of order 1
allocate( mm(mknot,0:ndeg,nage) )
mm = 0.d0
do i = 1,nage
   iage = iiage(i)
   do j = 1,mknot-1
      if(tvec(j) == tvec(j+1))cycle
      if(iage.ge.tvec(j).and.iage.lt.tvec(j+1)) mm(j,0,i)=1.d0
   end do
end do
! higher order basis functions
do ii = 1, ndeg
   do ia = 1,nage
      age = iiage(ia)
      do j = 1,mknot-ii
      if( tvec(j+ii) > tvec(j) )      mm(j,ii,ia) =                        &
&              ( age - tvec(j) ) / ( tvec(j+ii)-tvec(j)) * mm(j,ii-1,ia)
      if(j+ii+1 > mknot) cycle
      if( tvec(j+ii+1) > tvec(j+1) )  mm(j,ii,ia) = mm(j,ii,ia) +          &
&     ( tvec(j+ii+1) - age) / ( tvec(j+ii+1) - tvec(j+1)) * mm(j+1,ii-1,ia)
      end do
   end do
end do

! write out file of covariables
bsplfile(9:9)=degree(ndeg)
if(nrr <10 )then
   write(bsplfile(11:11),'(i1)') nrr
else
   write(bsplfile(10:11),'(i2)') nrr
end if
open(11,file=bsplfile,status='unknown',form='formatted')
write( fmt(2:3),'(i2)' ) nrr
do ia = 1, nage
   write(11,fmt) mm(:nrr,ndeg,ia), iiage(ia)
end do
close(11)

! write out a summary
  write(*,*)'number of ages            = ',nage
  write(*,*)' ... range                  ',iiage(1),iiage(nage)
  write(*,*)'degree of polynomial      = ',ndeg
  write(*,*)'number of knots           = ',nknot
  write(*,*)'number of RR coefficients = ',nrr
  if( equal ) write(*,*)'Uniform spline with auxiliary knots outside range'
  write(*,'(1x,a,1h",a,1h")')'output file               : ',bsplfile

  contains

! ==================
  subroutine cmdline
! ==================

! routine to evaluate command line options
  integer           :: iargc, i, k, l
  character(len=30) :: arg
  external getarg

! default value for options
  verbose =.false.; degset=.false.; kntset=.false.; equal=.false.; ages =.false.

  k = iargc() ! no. of arguments
  do i = 1, k
     call getarg (i,arg);  arg = adjustl(arg)
     if(arg(1:2) == '-v' .or. arg(1:2) == '-V') then
        verbose = .true.
     else if(arg(1:3) == '-eq' .or. arg(1:3) == '-EQ' ) then
        equal=.true.; read(arg(4:),*)nknot;  kntset=.true.
     else if(arg(1:2) == '-a' .or. arg(1:2) == '-A')then
        ages = .true.
     else if( arg(1:2) == '-d' .or. arg(1:2) == '-D') then
        read(arg(3:),*) ndeg; if(ndeg <1 .or. ndeg >3 ) stop 'wrong degree'
        degset = .true.
     else if( arg(1:2) =='-h' .or. arg(1:2) == '-H') then
        write(*,*)'Valid command line options :'
        write(*,'(a,t10,a)')'-h','produces this message'
        write(*,'(a,t10,a)')'-u','how to use it'
        write(*,'(a,t10,a)')'-v','verbose screen output'
        write(*,'(a,t10,a)')'-a','specify range of ages & step size'
        write(*,'(a,t10,a)')'-dx','for x=1,2,3 : degree of spline'
        write(*,'(a,t10,a)')'-eqx','for x=1,2,... : uniform spline with x knots'
     else if( arg(1:2) == '-u' .or. arg(1:2) == '-U' ) then
        write(*,*)'Program BSPLINES : set up B-spline functions    '
        write(*,*)'------------------------------------------------'
        write(*,*)'Input: file "ages.out"'
        write(*,*)'Input: File "knots" with user-defined knots (one per line)'
        write(*,*)'Input: Options via command line or interactively'
        write(*,*)'        - degree of spline (1=linear, 2=quadratic, 3=cubic)'
        write(*,*)'        - absence of "knots" selects a uniform spline with'
        write(*,*)'          auxiliary knots equi-distant outside range'
        write(*,*)'        - for "ages.out" absent, user is prompted for age'
        write(*,*)'          range and step size ';  stop
     end if
  end do
  if( verbose ) then
     write(*,*)' **** Command line options ***'
     if(degset) write(*,'(1x,a,t30,1h=,i12)')'Degree of spline',ndeg
     if(equal) write(*,'(1x,a,t30,1h=,i12,1x,a)')'Uniform spline with', &
&                                                      nknot,'knots'
  end if
  return
  end subroutine cmdline

  end program Bsplines