! --------------------------------------------------------------------------------
! Program FitPlane
! --------------------------------------------------------------------------------
!  Usage:  < xyz.H  FitPlane.x m1= m2= o1_grid= o2_grid= 
!                              d1_grid= d2_grid=            > fullplane.H
! --------------------------------------------------------------------------------
!  Description: From a collection of (x,y,z) data, computes the least squares
!               best-fit plane, given by z = Ax + By + C.  Outputs the values of 
!               that plane evaluated over a user defined output grid.
! --------------------------------------------------------------------------------
!  Input:  xyz.H - a SEPLib history file containing nd (x,y,z) triplets.
! --------------------------------------------------------------------------------
!  Output:  fullplane.H - The least squares best-fit plane to the (x,y,z) data,
!                         evaluated on a regular grid (user defined).
! --------------------------------------------------------------------------------
! Parameters:  m1 - Output model grid size, 1-axis
!              m2 - Output model grid size, 2-axis
!              o1_grid - Output model origin, 1-axis
!              o2_grid - Output model origin, 2-axis
!              d1_grid - Output model grid spacing, 1-axis
!              d2_grid - Output model grid spacing, 2-axis
!              d2_grid - Output model grid spacing, 2-axis
! --------------------------------------------------------------------------------

program FitPlane
  use sep
  use planefit_mod
  use cgstep_mod
  use dottest
  use solver_smp_mod

  implicit none

  integer                         :: nd               ! Number of (x,y,z) triplets
  integer                         :: m1, m2 
  integer                         :: i1, i2
  real                            :: om1, om2, dm1, dm2
  real                            :: x, y
  real                            :: dot1(2), dot2(2)         ! Dot product test
  real, dimension (:), pointer    :: coeff                    ! The model: [A B C]'
  real, dimension (:,:), pointer  :: xyz                      ! Input data
  real, dimension (:,:), pointer  :: planefull                ! Output grid

  call sep_init (SOURCE)

  ! ------------------------------------------------------------------------------
  !  Parameter reading
  ! ------------------------------------------------------------------------------
  call from_history( "n1", nd)
  call from_param( "m1",      m1);   call from_param( "m2",      m2)
  call from_param( "o1_grid", om1);  call from_param( "o2_grid", om2)
  call from_param( "d1_grid", dm1);  call from_param( "d2_grid", dm2)

  allocate( xyz(nd,3), planefull(m1,m2), coeff(3) )

  call sep_read( xyz )

  call planefit_init( xyz(:,1:2), nd)

  ! ------------------------------------------------------------------------------
  ! !!!!!!!!!!!!!!!!!!!!!! HEY, DIRECT YOUR GAZE HERE!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  ! ------------------------------------------------------------------------------
  !  USE solver_smp() TO ESTIMATE THE PLANE COEFFICIENTS.
  ! ------------------------------------------------------------------------------
!  !  call solver_smp(m=?,d=?,Fop=planefit_lop,stepper=?,niter=?, verb=.true.)


  ! ------------------------------------------------------------------------------
  !  Run Dot product test and write estimated coefficients to the screen.
  ! ------------------------------------------------------------------------------
  dot1=0; dot2=0
  write(0,*) '--- PLANEFIT DOT PRODUCT TEST -------------------------------------'
  call dot_test( planefit_lop, 3, nd, dot1, dot2)
  write(0,*) '-------------------------------------------------------------------'
  write(0,*) '--- THESE ARE YOUR PLANE COEFFICIENTS -----------------------------'
  write(0,*) coeff
  write(0,*) '-------------------------------------------------------------------'

  ! ------------------------------------------------------------------------------
  ! !!!!!!!!!!!!!!!!!!!!!! HEY, DIRECT YOUR GAZE HERE!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  ! ------------------------------------------------------------------------------
  !  USING THE PLANE COEFFICIENTS AND THE GRID PARAMETERS, COMPUTE THE VALUE OF
  !  THE PLANE OVER THE FULL GRID AND STORE IN PLANEFULL. 
  ! ------------------------------------------------------------------------------


  ! ------------------------------------------------------------------------------

  call to_history( "n1", m1);   call to_history( "n2", m2)
  call to_history( "o1", om1);  call to_history( "o2", om2)
  call to_history( "d1", dm1);  call to_history( "d2", dm2)

  call sep_write( planefull )

  call sep_close()
  deallocate( planefull, xyz, coeff )
!!! TA Solution
  call planefit_close() ! Just to free memory
!!!
  call exit( 0)

end program FitPlane
