c23456789012345678901234567890123456789012345678901234567890123456789012
c
c program to prepare data files for use by splining routine
c***********************************************************************
c
      implicit double precision(a-h,o-z)
      real xchi,plev,df
      dimension ahatm(3,4),qhat(4),qbing(4,4),sigma(3,3)
      character filnam*100
1000  format(a)
c
      write (*,*) 'Output file name? '
	  read (*,1000) filnam
	  open (unit=12,file=filnam)
c this code would logically apply if hrcl8.for did not hard code 7.815 at line 224
c	  write (*,*) 'Which method will you use to choose lambda?'
c	  write (*,*) '0. Specify a specific lambda.'
c	  write (*,*) '1. Cross validation.'
c	  write (*,*) '2. Pass through the known confidence regions.'
c	  write (*,*) '3. Discrepancy.'
c	  write (*,*) 'Enter 0, 1, 2, or 3. '
c	  read (*,*) ilam
c	  if (ilam.eq.2) then
c	    write (*,*) 'Spline together what level confidence regions (e.g. .95)? '
c	    read (*,*) plev
c      endif
      plev=.95
      write (*,*) 'How many rotation files? '
	  read (*,*) nrot
	  write (12,*) nrot
c
      do 1 irot = 1, nrot	  
      write (*,*) 'Enter rotation file name for rotation number ',irot
      read (*,1000) filnam
      open (unit=11,file=filnam,status='old',err=1)
	  write (*,*) '  Age for rotation number ',irot
	  read (*,*) age
c
c read 2 first comment lines
      read(11,1000) filnam
      read(11,1000) filnam
c read bestfit rotation: alat,alon and angle
      read(11,1000) filnam
      read(11,*) alat,along,rho
      write(12,*) age,alat,along,rho
c=====translate best lat & lon & angle into a quaternion
      call trans3(alat,along,rho,qhat)
c read confidence level
      read(11,1000) filnam
      read(11,*) dum
c read hatkap and degree of freedom
      read(11,1000) filnam
      read(11,*) hatkap,df
      if (df.eq.0.) df=10000.
      if (hatkap.eq.0.) hatkap=1
c read nbr of points, sections
      read(11,1000) filnam
      read(11,*) dum
c read rotation matrix
      read(11,1000) filnam
	  read(11,*) dum
	  read(11,*) dum
	  read(11,*) dum
c read covariance matrix
      read(11,1000) filnam
	  read(11,*) dum
	  read(11,*) dum
	  read(11,*) dum
c read inverse of covariance matrix - unscaled: H11.2
      read(11,1000) filnam
      do 15 i=1,3
   15 read(11,*) (sigma(i,j),j=1,3)
c
      write (*,*) '  Do you want kappahat to be set to 1?'
      write (*,*) '          1=no    2=yes'
      read (*,*) ichoice
      if (ichoice.eq.2) hatkap=1
	  xchi=7.815
c      xchi=1
c	  if (ilam.eq.2) then
c        if (df.ge.1000.or.ichoice.eq.2) then
c          call xidch(plev,3.,xchi,ier)
c          if (ier.ne.0) write(6,*) 'subroutine xidch--ier ',ier
c        else
c          call xidf(plev,3.,df,xchi,ier)
c          if (ier.ne.0) write(6,*) 'subroutine xidf--ier ',ier
c          xchi=3.*xchi
c        endif
c      endif
      xchi=xchi/hatkap
c
      ahatm(1,1)=-qhat(2)
      ahatm(1,2)=qhat(1)
      ahatm(1,3)=qhat(4)
      ahatm(1,4)=-qhat(3)
      ahatm(2,1)=-qhat(3)
      ahatm(2,2)=-qhat(4)
      ahatm(2,3)=qhat(1)
      ahatm(2,4)=qhat(2)
      ahatm(3,1)=-qhat(4)
      ahatm(3,2)=qhat(3)
      ahatm(3,3)=-qhat(2)
      ahatm(3,4)=qhat(1)
      do 20 i=1,4
      do 20 j=1,4
      qbing(i,j)=0.
      do 20 k1=1,3
      do 20 k2=1,3
 20   qbing(i,j)=qbing(i,j)+4.*ahatm(k1,i)*sigma(k1,k2)*ahatm(k2,j)/xchi
c
      do 21 i=1,4
 21   write(12,*) (qbing(i,j),j=1,4)
 1    continue
      stop
      end