/*

# Copyright (c) 2005, 2008 Andrew Fraser
# This file is part of HMM_DS_Code.
#
# HMM_DS_Code is free software; you can redistribute it and/or modify it under
# the terms of the GNU General Public License as published by the Free Software
# Foundation; either version 2 of the License, or (at your option) any later
# version.
#
# HMM_DS_Code is distributed in the hope that it will be useful, but
# WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
# General Public License for more details.
#
# You should have received a copy of the GNU General Public License along with
# this program; if not, write to the Free Software Foundation, Inc.,
# 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
*/

#include <stdio.h>
#include <math.h>

struct ParPac {
  double s,b,r;
};

extern int Lsteps(double *y,
		   struct ParPac *pac,
		   double ts,
		   int Nt,
		   float result[][3]);
extern int Ltansteps(
    double *y,
    struct ParPac *pac,
    double ts,
    int Nt,
    double result[][12]);

int main(int argc, char *argv[]){

  double s = 1.26;
  double b = 0.08;
  double r = 66.9;
  double ts = 7.5;
  int Nt = 3;
  struct ParPac pac;
  double ic[3] = { 1,    1,   1};

  
  /* Pack the parameter pac */
  pac.s =s;  pac.b=b;  pac.r=r;

  /* Test Lsteps */
  {
    int i;
    float result[Nt][3];       /* Avoid simple array to make python */
    Lsteps(ic,&pac,ts,Nt,result);
    fprintf(stderr,"After Lsteps() in Ltest.c. Nt=%d\n",Nt);
    for(i=0;i<Nt;i++){
      fprintf(stderr,"%f %f  %f\n",result[i][0],result[i][1],result[i][2]);
    }
    for(i=0;i<3;i++)
      ic[i] = result[Nt-1][i];
  }
  /* Test Ltansteps */
  {
    int i,j;
    double y[12] = {ic[0],ic[1],ic[2],1,0,0,0,1,0,0,0,1};
    double result[Nt][12];
    ts = 0.15;
    Ltansteps(y,&pac,ts,Nt,result);
    fprintf(stderr,"Results of Ltansteps() Ltest.c \n");
    for(i=0;i<Nt;i++){
      printf("\n%7.3lf %7.3lf  %7.3lf\n\n",result[i][0],
	     result[i][1],result[i][2]);
      for(j=1;j<4;j++){
	printf("          %7.3lf %7.3lf  %7.3lf\n",
	       result[i][3*j],result[i][3*j+1],result[i][3*j+2]);
      }
    }
    }
  return 0;
}