diff -c -r -b -B cosmomc/params.ini cosmomc_sdss/params.ini
*** cosmomc/params.ini	2004-09-29 21:09:54.000000000 +0100
--- cosmomc_sdss/params.ini	2004-10-28 22:32:30.642309947 +0100
***************
*** 31,36 ****
--- 31,38 ----
  use_CMB = T
  use_HST = T
  use_mpk = F
+ use_2dF = F
+ use_sdss = F
  use_clusters = F
  use_BBN = F
  use_Age_Tophat_Prior = T
diff -c -r -b -B cosmomc/source/calclike.f90 cosmomc_sdss/source/calclike.f90
*** cosmomc/source/calclike.f90	2004-09-29 21:05:38.000000000 +0100
--- cosmomc_sdss/source/calclike.f90	2004-10-28 22:48:03.341921967 +0100
***************
*** 2,8 ****
   use CMB_Cls
   use cmbtypes
   use cmbdata
!  use mpk
   use Random
   use settings
   use ParamDef
--- 2,9 ----
   use CMB_Cls
   use cmbtypes
   use cmbdata
! ! use mpk
!  use twodf
   use Random
   use settings
   use ParamDef
diff -c -r -b -B cosmomc/source/cmbtypes.f90 cosmomc_sdss/source/cmbtypes.f90
*** cosmomc/source/cmbtypes.f90	2004-09-29 21:05:38.000000000 +0100
--- cosmomc_sdss/source/cmbtypes.f90	2004-10-28 22:03:17.061337253 +0100
***************
*** 17,25 ****
    integer, parameter :: num_matter_power = 49 !number of points computed in matter power spectrum
    real, parameter    :: matter_power_minkh =  3.16227766E-03 !minimum value of k/h to store
    real, parameter    :: matter_power_dlnkh = 0.143911568     !log spacing in k/h
!   real, parameter    :: matter_power_maxz = 0.    !6.0
!   integer, parameter :: matter_power_lnzsteps = 1 !20
! 
  
          !Note these are the interpolated/extrapolated values. The k at which matter power is computed up to 
          !by CAMB is set in CMB_Cls_xxx with, e.g. P%Transfer%kmax = 0.6 (which is enough for 2dF)
--- 17,28 ----
    integer, parameter :: num_matter_power = 49 !number of points computed in matter power spectrum
    real, parameter    :: matter_power_minkh =  3.16227766E-03 !minimum value of k/h to store
    real, parameter    :: matter_power_dlnkh = 0.143911568     !log spacing in k/h
! !sdss
! !  real, parameter    :: matter_power_maxz = 0.    !6.0
! !  integer, parameter :: matter_power_lnzsteps = 1 !20
!   real, parameter    :: matter_power_maxz = 0.17    !6.0
!   integer, parameter :: matter_power_lnzsteps = 2 !20
! !end sdss
  
          !Note these are the interpolated/extrapolated values. The k at which matter power is computed up to 
          !by CAMB is set in CMB_Cls_xxx with, e.g. P%Transfer%kmax = 0.6 (which is enough for 2dF)
***************
*** 323,335 ****
       integer i
     
       x = log(kh/matter_power_minkh) / matter_power_dlnkh
       if (x < 0 .or. x >= num_matter_power-1) then
!         write (*,*) ' k/h out of bounds in MatterPowerAt (',kh,')'
!         stop
!      end if
       i = int(x)
       d = x - i
       MatterPowerAt = exp(log(T%matter_power(i+1,1))*(1-d) + log(T%matter_power(i+2,1))*d)
       !Just do linear interpolation in logs for now..
       !(since we already cublic-spline interpolated to get the stored values)
  
--- 326,344 ----
       integer i
     
       x = log(kh/matter_power_minkh) / matter_power_dlnkh
+ 
+ ! SDSS window function need a wide range. Setting nul power spectrum
+ !from extremal values does not change the chi2 significantly.
+ 
       if (x < 0 .or. x >= num_matter_power-1) then
! !        write (*,*) ' k/h out of bounds in MatterPowerAt (',kh,')'
! !        stop
!         MatterPowerAt=0.
!      else
          i = int(x)
          d = x - i
          MatterPowerAt = exp(log(T%matter_power(i+1,1))*(1-d) + log(T%matter_power(i+2,1))*d)
+      endif
       !Just do linear interpolation in logs for now..
       !(since we already cublic-spline interpolated to get the stored values)
  
***************
*** 355,364 ****
          stop
       end if
       x = log(kh/matter_power_minkh) / matter_power_dlnkh
       if (x < 0 .or. x >= num_matter_power-1) then
!         write (*,*) ' k/h out of bounds in MatterPowerAt_Z (',kh,')'
!         stop
!      end if
  
       iz = int(y)
       dz = y - iz
--- 362,377 ----
          stop
       end if
       x = log(kh/matter_power_minkh) / matter_power_dlnkh
+ 
+ ! SDSS window function need a wide range. Setting nul power spectrum
+ !from extremal values does not change the chi2 significantly.
+ ! if (dz...) to respect array bound
+ 
       if (x < 0 .or. x >= num_matter_power-1) then
! !        write (*,*) ' k/h out of bounds in MatterPowerAt_Z (',kh,')'
! !        stop
!         MatterPowerAt_Z=0.
!      else
  
          iz = int(y)
          dz = y - iz
***************
*** 366,376 ****
       i = int(x)
       d = x - i
  
!      mup = log(T%matter_power(i+1,iz+2))*(1-d) + log(T%matter_power(i+2,iz+2))*d
!      mdn = log(T%matter_power(i+1,iz+1))*(1-d) + log(T%matter_power(i+2,iz+1))*d
  
       MatterPowerAt_Z = exp(mdn*(1-dz) + mup*dz)
  
     end function MatterPowerAt_Z
  
  
--- 379,393 ----
          i = int(x)
          d = x - i
  
!         mup=0.
!         mdn=0.
!         if (dz.ne.0.) mup = log(T%matter_power(i+1,iz+2))*(1-d) + log(T%matter_power(i+2,iz+2))*d
!         if (dz.ne.1.) mdn = log(T%matter_power(i+1,iz+1))*(1-d) + log(T%matter_power(i+2,iz+1))*d
  
          MatterPowerAt_Z = exp(mdn*(1-dz) + mup*dz)
  
+      endif
+ 
     end function MatterPowerAt_Z
  
  
diff -c -r -b -B cosmomc/source/driver.F90 cosmomc_sdss/source/driver.F90
*** cosmomc/source/driver.F90	2004-10-14 14:36:03.000000000 +0100
--- cosmomc_sdss/source/driver.F90	2004-10-28 22:57:00.938273081 +0100
***************
*** 26,31 ****
--- 26,33 ----
          integer status          
          real delta_loglike
  
+ 	logical :: Use_mpk
+ 
  #ifdef MPI
          double precision intime
          integer ierror
***************
*** 141,149 ****
          if (Ini_Read_Logical('cmb_hyperparameters', .false.)) &
              call DoStop( 'Hyperparameters not supported any more')
  
-         if (Ini_Read_String('use_2dF') /= '') stop 'use_2dF now replaced with use_mpk'
          Use_Clusters = Ini_Read_Logical('use_clusters',.false.)
!         Use_mpk = Ini_Read_Logical('use_mpk',.false.) ! matter power spectrum, incl 2dF
          Use_HST = Ini_Read_Logical('use_HST',.true.)
          Use_BBN = Ini_Read_Logical('use_BBN',.false.)
          Use_Age_Tophat_Prior= Ini_Read_Logical('use_Age_Tophat_Prior',.true.)
--- 143,160 ----
          if (Ini_Read_Logical('cmb_hyperparameters', .false.)) &
              call DoStop( 'Hyperparameters not supported any more')
  
          Use_Clusters = Ini_Read_Logical('use_clusters',.false.)
!         
! 	Use_2dF = Ini_Read_Logical('use_2dF',.false.)
! 	Use_sdss = Ini_Read_Logical('use_sdss',.false.)
! 	write(*,*)' ***********************************************************'
! 	write(*,*) 'SDSS is now implemented in COSMOMC!'
! 	write(*,*) 'this patch is for testing purpose and you should not use it!'
! 	write(*,*)' ***********************************************************'
! 	read(*,*)
! 	Use_mpk=.false.
!        
! 
          Use_HST = Ini_Read_Logical('use_HST',.true.)
          Use_BBN = Ini_Read_Logical('use_BBN',.false.)
          Use_Age_Tophat_Prior= Ini_Read_Logical('use_Age_Tophat_Prior',.true.)
***************
*** 153,159 ****
          Use_min_zre = Ini_Read_Double('use_min_zre',0.d0) 
  
          ! use_LSS = Use_2dF .or. Use_Clusters .or. Use_WeakLen
!         use_LSS = Use_mpk .or. Use_Clusters .or. Use_WeakLen
  
          Temperature = Ini_Read_Real('temperature',1.)
          
--- 164,170 ----
          Use_min_zre = Ini_Read_Double('use_min_zre',0.d0) 
  
          ! use_LSS = Use_2dF .or. Use_Clusters .or. Use_WeakLen
!         use_LSS = Use_mpk .or. Use_Clusters .or. Use_WeakLen .or. Use_2dF .or. Use_sdss
  
          Temperature = Ini_Read_Real('temperature',1.)
          
diff -c -r -b -B cosmomc/source/Makefile cosmomc_sdss/source/Makefile
*** cosmomc/source/Makefile	2004-09-30 20:07:02.000000000 +0100
--- cosmomc_sdss/source/Makefile	2004-10-28 22:52:57.790979291 +0100
***************
*** 94,100 ****
  DISTFILES = utils.o settings.o GetDist.o
  
  OBJFILES= WMAP_likelihood.o\
! 	utils.o settings.o cmbtypes.o cmbdata.o WeakLen.o twodf.o mpk.o supernovae.o\
  	$(CLSFILE) paramdef.o $(PROPOSE) $(PARAMETERIZATION) calclike.o \
  	conjgrad_wrapper.o EstCovmat.o postprocess.o MCMC.o driver.o
  
--- 94,100 ----
  DISTFILES = utils.o settings.o GetDist.o
  
  OBJFILES= WMAP_likelihood.o\
! 	utils.o settings.o cmbtypes.o cmbdata.o WeakLen.o twodf.o supernovae.o\
  	$(CLSFILE) paramdef.o $(PROPOSE) $(PARAMETERIZATION) calclike.o \
  	conjgrad_wrapper.o EstCovmat.o postprocess.o MCMC.o driver.o
  
diff -c -r -b -B cosmomc/source/Makefile_intel cosmomc_sdss/source/Makefile_intel
*** cosmomc/source/Makefile_intel	2004-10-13 13:50:20.000000000 +0100
--- cosmomc_sdss/source/Makefile_intel	2004-10-28 22:53:59.200343729 +0100
***************
*** 46,52 ****
  DISTFILES = ../camb/inifile.o utils.o settings.o GetDist.o
  
  OBJFILES= $(CAMBFILES) WMAP_likelihood.o\
! 	utils.o settings.o cmbtypes.o cmbdata.o WeakLen.o twodf.o mpk.o supernovae.o\
  	$(CLSFILE) paramdef.o $(PROPOSE) $(PARAMETERIZATION) calclike.o \
  	conjgrad_wrapper.o EstCovmat.o postprocess.o MCMC.o driver.o
  
--- 46,52 ----
  DISTFILES = ../camb/inifile.o utils.o settings.o GetDist.o
  
  OBJFILES= $(CAMBFILES) WMAP_likelihood.o\
! 	utils.o settings.o cmbtypes.o cmbdata.o WeakLen.o twodf.o supernovae.o\
  	$(CLSFILE) paramdef.o $(PROPOSE) $(PARAMETERIZATION) calclike.o \
  	conjgrad_wrapper.o EstCovmat.o postprocess.o MCMC.o driver.o
  
diff -c -r -b -B cosmomc/source/twodf.f90 cosmomc_sdss/source/twodf.f90
*** cosmomc/source/twodf.f90	2004-09-29 21:05:38.000000000 +0100
--- cosmomc_sdss/source/twodf.f90	2004-02-12 15:06:45.000000000 +0000
***************
*** 1,5 ****
! !Module for using 147k 2dF data 
  !See http://www.roe.ac.uk/~wjp/2dFGRS/2dF_Pk.html - adapted for f90
  
  module twodf
  use settings
--- 1,9 ----
! !Module for using 147k 2dF data and SDSS data
  !See http://www.roe.ac.uk/~wjp/2dFGRS/2dF_Pk.html - adapted for f90
+ !See http://www.hep.upenn.edu/~max/sdss.html - adapted for f90
+ 
+ ! version february 2004.
+ ! TODO: introduce bias as fast parameter; implement corrections from HaloFit.
  
  module twodf
  use settings
***************
*** 8,17 ****
--- 12,54 ----
  
   !Note all units are in k/h here
  
+  logical :: use_2dF = .false.
+ 
   integer, parameter :: twodf_points = 32
   real :: twodf_P(twodf_points), twodf_k(twodf_points), twodf_Ninv(twodf_points,twodf_points)
   logical :: do_twodf_init  = .true.
  
+  real, parameter :: twodf_redshift = 0.17
+ 
+ !SDSS declaration
+ 
+  logical :: use_sdss = .false.
+ 
+ !To use hyper parameters between twodf and sdss
+  logical :: use_lss_hyper=.false.
+ 
+  integer, parameter :: sdss_bands_full = 22
+  integer, parameter :: sdss_kbands_full = 97
+  
+ !The window functions are defined up to kh=100 for kbands=97. The
+ !observed power spectrum is measured up to kh=0.3 for bands=22. Linear
+ !regime is trusted up to k/h=0.15 or bands=17 (k/h=0.2 for 19). The
+ !extrapolated power spectrum is arbitrarily set to zero for k/h>3 (see
+ !cmbtype.f90)
+ 
+  integer, parameter :: sdss_bands = 17 
+  integer, parameter :: sdss_kbands = 97
+  real, parameter :: sdss_redshift = 0.1
+  real, dimension(sdss_bands) :: sdss_P
+  real, dimension(sdss_kbands) :: sdss_k
+  real, dimension(sdss_bands,sdss_kbands) :: sdss_W
+  logical :: do_sdss_init  = .true.
+ 
+  !Ln(Gamma(ndata/2 + 1)) for hyperparameters
+  real :: Twodfo2p1,Sdsso2p1
+  real :: LnGammaTwodf,LnGammaSdss
+  
+ 
  contains
  
   subroutine twodf_init
***************
*** 26,34 ****
      end do
      close(tmp_file_unit)
      call ReadMatrix('data/twodf_june02_matrix.dat',twodf_Ninv,twodf_points,twodf_points)
! 
     do_twodf_init = .false.
! 
      read (tmp_file_unit,*, end = 300, err=300) model 
  
  200 stop 'Error reading 2dF file'
--- 63,72 ----
      end do
      close(tmp_file_unit)
      call ReadMatrix('data/twodf_june02_matrix.dat',twodf_Ninv,twodf_points,twodf_points)
!     Twodfo2p1 = real(twodf_points)/2. + 1.
!     LnGammaTwodf = gammln(Twodfo2p1)
      do_twodf_init = .false.
!     if (Feedback > 2) write(*,*) 'Ln[Gamma(1+twodf_points/2)]= ',LnGammaTwodf
      read (tmp_file_unit,*, end = 300, err=300) model 
  
  200 stop 'Error reading 2dF file'
***************
*** 74,80 ****
        do imu = 1, nint2
           rmu = -1 + (imu-0.5)*dmu
           yk = sqrt(rk2 + xk2 - 2*xk*rk*rmu)
!          sp = sp + MatterPowerAt(Theory, yk)
        end do
        s = s + sp*xk2*twodf_Wksq(xk)
        xk = rkmax; xk2= xk**2
--- 112,119 ----
        do imu = 1, nint2
           rmu = -1 + (imu-0.5)*dmu
           yk = sqrt(rk2 + xk2 - 2*xk*rk*rmu)
! !         sp = sp + MatterPowerAt(Theory, yk)
!          sp = sp + MatterPowerAt_Z(Theory, yk,twodf_redshift)
        end do
        s = s + sp*xk2*twodf_Wksq(xk)
        xk = rkmax; xk2= xk**2
***************
*** 97,103 ****
           do imu = 1, nint2
              rmu = -1 + (imu-0.5)*dmu
              yk = sqrt(rk2+xk2-2*xk*rk*rmu)
!             sp = sp + MatterPowerAt(Theory, yk)
           end do
           s = s+ sp*xk2*twodf_Wksq(xk)
        end do
--- 136,143 ----
           do imu = 1, nint2
              rmu = -1 + (imu-0.5)*dmu
              yk = sqrt(rk2+xk2-2*xk*rk*rmu)
! !            sp = sp + MatterPowerAt(Theory, yk)
!             sp = sp + MatterPowerAt_Z(Theory, yk,twodf_redshift)
           end do
           s = s+ sp*xk2*twodf_Wksq(xk)
        end do
***************
*** 113,119 ****
   function LSS_2dflike(CMB, Theory)
     Type (CMBParams) CMB
     Type (CosmoTheory) Theory
!    real LSS_2dflike
     real th(twodf_points), sum1, chisq
     real Ninvth(twodf_points), Ninvd(twodf_points)
     integer i
--- 153,159 ----
   function LSS_2dflike(CMB, Theory)
     Type (CMBParams) CMB
     Type (CosmoTheory) Theory
!    real, dimension(2) :: LSS_2dflike
     real th(twodf_points), sum1, chisq
     real Ninvth(twodf_points), Ninvd(twodf_points)
     integer i
***************
*** 139,148 ****
      chisq = sum(twodf_P*Ninvd) - sum(th*Ninvd)**2/sum1 + log(sum1)
     
  
!     LSS_2dflike = chisq/2
  
      if (Feedback>1) write(*,*) '2df chi-sq: ', chisq
     
   end function LSS_2dflike
  
  end module twodf
--- 179,336 ----
      chisq = sum(twodf_P*Ninvd) - sum(th*Ninvd)**2/sum1 + log(sum1)
     
  
!     LSS_2dflike(1) = chisq/2.
!     LSS_2dflike(2) = chisq/2. - log(sum1)/2.
  
      if (Feedback>1) write(*,*) '2df chi-sq: ', chisq
     
   end function LSS_2dflike
  
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! The SDSS stuff
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ 
+  subroutine sdss_init   
+           
+    integer i,j,iopb
+    real keff,klo,khi,beff
+    real, dimension(sdss_bands_full,sdss_kbands_full) :: sdss_Wfull
+    real, dimension(sdss_kbands_full) :: sdss_kfull
+    real, dimension(sdss_bands) :: sdss_sdev
+    character(80) :: dummychar
+    
+    call ReadVector('data/sdss_kbands.txt',sdss_kfull,sdss_kbands_full)
+    sdss_k(1:sdss_kbands)=sdss_kfull(1:sdss_kbands)
+ 
+    if (Feedback > 0) then 
+       write(*,*) 'reading: sdss data'
+       write(*,*) 'kbands Windows truncated at k/h=  ',sdss_k(sdss_kbands)      
+    endif
+    call OpenTxtFile('data/sdss_measurements.txt', tmp_file_unit)
+    sdss_P=0.
+    read (tmp_file_unit,*)dummychar
+    read (tmp_file_unit,*)dummychar
+    do i =1, sdss_bands
+       read (tmp_file_unit,*, iostat=iopb) keff,klo,khi,sdss_P(i),sdss_sdev(i),beff
+    end do
+    close(tmp_file_unit)
+ 
+    if (Feedback > 0) write(*,*) 'bands truncated at keff=  ',keff
+    
+    call ReadMatrix('data/sdss_windows.txt',sdss_Wfull,sdss_bands_full,sdss_kbands_full)
+    sdss_W(1:sdss_bands,1:sdss_kbands)=sdss_Wfull(1:sdss_bands,1:sdss_kbands)
+    
+   
+    ! Prewhiten so that WLOG N=I:
+    do i=1,sdss_bands
+       sdss_P(i) = sdss_P(i)/sdss_sdev(i)
+       do j=1,sdss_kbands
+          sdss_W(i,j) = sdss_W(i,j)/sdss_sdev(i)
+       end do
+    end do
+    ! Now chi2 = (obs-W p)^t (obs-W p) = |obs-Wp|^2.
+    Sdsso2p1 = real(sdss_bands)/2. + 1.
+    LnGammaSdss = gammln(Sdsso2p1)
+    do_sdss_init = .false.
+    if (Feedback > 2) write(*,*) 'Ln[Gamma(1+sdss_bands/2)]= ',LnGammaSdss
+    if (iopb.ne.0) then
+       stop 'Error reading sdss file'
+    endif
+  
+  end subroutine sdss_init
+ 
+  
+  function LSS_sdsslike(CMB, Theory)
+   Type (CMBParams) CMB
+   Type (CosmoTheory) Theory
+   real, dimension(sdss_kbands) :: sdss_Pth
+   real, dimension(sdss_bands) :: sdss_WPth
+   real :: chisq,normV
+   real, dimension(2) :: LSS_sdsslike
+   integer :: i
+     
+   if (do_sdss_init) call sdss_init
+ 
+   do i=1, sdss_kbands
+      sdss_Pth(i)=MatterPowerAt_Z(Theory,sdss_k(i),sdss_redshift)     
+   end do
+   
+   sdss_WPth = matmul(sdss_W,sdss_Pth)
+   
+   !with analytic marginalization over normalization nuisance (flat prior)
+   normV = sum(sdss_WPth*sdss_WPth)
+   chisq = sum(sdss_P*sdss_P) - sum(sdss_WPth*sdss_P)**2/normV + log(normV)
+   
+   !without analytic marginalization
+   !  chisq = sum((sdss_P(:) - sdss_WPth(:))**2)
+   
+   LSS_sdsslike(1) = chisq/2.
+   LSS_sdsslike(2) = chisq/2. - log(normV)/2.
+ 
+   if (Feedback>1) write(*,*) 'sdss chi-sq:', chisq
+    
+ end function LSS_sdsslike
+ 
+ 
+ function gammln(xx)
+   ! from numrec
+   REAL gammln,xx
+   INTEGER j
+   DOUBLE PRECISION ser,stp,tmp,x,y,cof(6)
+   SAVE cof,stp
+   DATA cof,stp/76.18009172947146d0,-86.50532032941677d0, &
+        24.01409824083091d0,-1.231739572450155d0,.1208650973866179d-2, &
+        -.5395239384953d-5,2.5066282746310005d0/
+   x=xx
+   y=x
+   tmp=x+5.5d0
+   tmp=(x+0.5d0)*log(tmp)-tmp
+   ser=1.000000000190015d0
+   do j=1,6
+      y=y+1.d0
+      ser=ser+cof(j)/y
+   enddo
+   gammln=tmp+log(stp*ser/x)
+   return
+ END function gammln
+ 
+ 
+ function LSSLnLike(CMB, Theory)
+   Type (CMBParams) CMB
+   Type (CosmoTheory) Theory
+   real LSSLnLike
+   real, dimension(2) :: twodflike,sdsslike
+   
+   LSSLnLike = 0.
+   
+   if (use_lss_hyper.and.use_2dF.and.use_sdss) then
+      twodflike = LSS_2dflike(CMB,Theory)
+      sdsslike = LSS_sdsslike(CMB,Theory)
+ 
+      LSSLnLike = Twodfo2p1*log(1.+twodflike(2)) + twodflike(1)-twodflike(2) - LnGammaTwodf &
+           + Sdsso2p1*log(1.+sdsslike(2)) + sdsslike(1)-sdsslike(2) - LnGammaSdss
+ 
+      if (Feedback>1) then
+         write(*,*) 'Effective hyperparam 2dF:  ',2.*real(twodf_points)/twodflike(2)
+         write(*,*) 'Effective hyperparam sdss: ',2.*real(sdss_bands)/sdsslike(2)        
+         write(*,*) 'LSS marginalized -ln(like):',LSSLnLike
+      endif
+ 
+      if (Feedback>2) then
+         write(*,*) 'Ln[Evidence ratio for hyperparams]: ',-LSSLnLike + twodflike(1) + sdsslike(1)
+      endif
+ 
+   else
+      if (use_2dF) then
+         twodflike = LSS_2dflike(CMB,Theory)
+         LSSLnLike = twodflike(1) + LSSLnLike
+      endif
+      if (use_sdss) then
+         sdsslike = LSS_sdsslike(CMB,Theory)
+         LSSLnLike = sdsslike(1) + LSSLnLike
+      endif
+   endif
+ 
+ end function LSSLnLike
+ 
  end module twodf
