 /*
  * Khoros: $Id: lmlde.c,v 1.2 1992/03/20 23:46:50 dkhoros Exp $
  */

#if !defined(lint) && !defined(SABER)
static char rcsid[] = "Khoros: $Id: lmlde.c,v 1.2 1992/03/20 23:46:50 dkhoros Exp $";
#endif

 /*
  * $Log: lmlde.c,v $
 * Revision 1.2  1992/03/20  23:46:50  dkhoros
 * VirtualPatch5
 *
  */

/*
 *----------------------------------------------------------------------
 *
 * Copyright 1992, University of New Mexico.  All rights reserved.
 * Permission to copy and modify this software and its documen-
 * tation only for internal use in your organization is hereby
 * granted, provided that this notice is retained thereon and
 * on all copies.  UNM makes no representations as to the sui-
 * tability and operability of this software for any purpose.
 * It is provided "as is" without express or implied warranty.
 * 
 * UNM DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE,
 * INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FIT-
 * NESS.  IN NO EVENT SHALL UNM BE LIABLE FOR ANY SPECIAL,
 * INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY OTHER DAMAGES WHAT-
 * SOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER
 * IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS
 * ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR PER-
 * FORMANCE OF THIS SOFTWARE.
 * 
 * No other rights, including, for example, the right to redis-
 * tribute this software and its documentation or the right to
 * prepare derivative works, are granted unless specifically
 * provided in a separate license agreement.
 *---------------------------------------------------------------------
 */

#include "unmcopyright.h"        /* Copyright 1992 by UNM */

/*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>  <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 >>>>
 >>>>         File Name: lmlde.c
 >>>>
 >>>>      Program Name: mlde
 >>>>
 >>>> Date Last Updated: Thu Mar  5 08:51:31 1992 
 >>>>
 >>>>          Routines: lmlde - the library call for mlde
 >>>>
 >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>   <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<*/


#include "vinclude.h"


/* -library_includes */
#define CLEANUP \
    {  /* cleanup workspace before leaving */ \
       if(y!=NULL)free(y); \
       if(dy!=NULL)free(dy); \
       if(dym!=NULL)free(dym); \
       if(dyt!=NULL)free(dyt); \
       if(yt!=NULL)free(yt); \
    } /* end of cleanup code */

typedef struct complex {
  double r,i;
} complex;
/* -library_includes_end */


/****************************************************************
*
* Routine Name: lmlde - library call for mlde
*
* Purpose:
*    
*    solves a system of linear ODE's
*    
*    

* Input:
*    
*    t0             initial time value for the solution.
*    
*    tf             final time value for the solution
*    
*    h              time step to be used to progress from t0 to tf.
*    
*    image          real or complex coefficient matrix.
*    
*    bimage         real or complex  vector  containing  the  constant
*                   driving functions of the system.
*    
*    ivimage        real or  complex  vector  containing  the  initial
*                   values of the system.
*    
*    

* Output:
*    
*    image          complex vector containing  solution  to  the  dif-
*                   ferential   equations.  Solution  is  interleaved.
*                   i.e., if matrix has four rows, then solution  will
*                   have  five elements per row.  The first element is
*                   the time, the remaining elements are the solutions
*                   to the system at various that time step.
*    
*    Return Value:  1 on success, 0 on failure.
*    
*    

*
* Written By: Jeremy Worley
*    
*    Jeremy Worley 05 Mar 1992 08:50 MST
*              Made certain implicitly declared functions explicit  in
*              both lmlde() and llde().
*    
*    

****************************************************************/


/* -library_def */
int lmlde(image,bimage,ivimage,t0,tf,h)
    struct xvimage *image,*bimage,*ivimage;
    float t0,tf,h;
/* -library_def_end */

/* -library_code */
{
  int numpts,i,j,rows1,cols1,rows2,cols2,type;
  double **mat,*m1,*m2,*m3;
  float *locs;
  char *program = "lmlde";
  struct xvimage *images[3];
  int proper_num_bands(), proper_num_images(), match_num_bands(), 
      lvcast(), llde();

  images[0] = image;
  images[1] = bimage;
  images[2] = ivimage;

  rows1 = (int)image->col_size;
  cols1 = (int)image->row_size;
  rows2 = (int)bimage->col_size;
  cols2 = (int)bimage->row_size;

/*
** check some stuff out
*/

  if(rows2!=cols1 || cols2!=1){
     fprintf(stderr,"%s:  improper dimension on b-vector.\n",program);
     return(0);
  }

  if(!proper_num_bands(program,image,1,FALSE)){
     fprintf(stderr,"%s:  image structure should have only one band.\n",
         program);
     return(0);
  }
 
  if(!proper_num_images(program,image,1,FALSE)){
     fprintf(stderr,"%s:  image structures should have only one image.\n",
         program);
     return(0);
  }

  if(!proper_num_images(program,bimage,1,FALSE)){
     fprintf(stderr,"%s:  image structures should have only one image.\n",
         program);
     return(0);
  }

  if(!proper_num_images(program,ivimage,1,FALSE)){
     fprintf(stderr,"%s:  image structures should have only one image.\n",
         program);
     return(0);
  }

  if(!match_num_bands(program,image,bimage,FALSE)){
     fprintf(stderr,"%s:  images structures must have same number of bands.\n",
         program);
     return(0);
  }

  if(!match_num_bands(program,image,ivimage,FALSE)){
     fprintf(stderr,"%s:  images structures must have same number of bands.\n",
         program);
     return(0);
  }

/*
** cast all of the data up to double or double complex
*/

   if(image->data_storage_type==VFF_TYP_COMPLEX || image->data_storage_type==
        VFF_TYP_DCOMPLEX || bimage->data_storage_type==VFF_TYP_COMPLEX ||
        bimage->data_storage_type==VFF_TYP_DCOMPLEX || 
        ivimage->data_storage_type==VFF_TYP_DCOMPLEX || 
        ivimage->data_storage_type==VFF_TYP_DCOMPLEX){
      if(!lvcast(images,(int)3,(int)VFF_TYP_DCOMPLEX)){
         fprintf(stderr,"%s:  Failure to cast up to double complex.\n",program);
         return(0);
      }
      type = 1;
   }else{
      if(!lvcast(images,(int)3,(int)VFF_TYP_DOUBLE)){
         fprintf(stderr,"%s:  Failure to cast up to double.\n",program);
         return(0);
      }
      type = 0;
   }

   numpts = (int)((tf-t0)/h) + 2;

/*
** process all of the data
*/

   m1 = (double *)(image->imagedata);
   m2 = (double *)(bimage->imagedata);
   m3 = (double *)(ivimage->imagedata);

   if((mat = (double **)malloc((cols1+1)*sizeof(double *)*(type+1)))
      ==NULL){
      fprintf(stderr,"%s:  No dice...malloc failed.\n",program);
      return(0);
   }

   if(!llde(mat,(double)t0,(double)tf,(double)h,m1,m2,m3,
      rows1,cols1,type,type,type)){
      fprintf(stderr,"%s:  Failed call to static routine.\n",program);
      free((char *)mat);
      return(0);
   }

   free((char *)image->imagedata);
   if((image->imagedata = (char *)malloc(numpts*cols1*sizeof(double)*(type+1)))
      ==NULL){
      fprintf(stderr,"%s:  No dice...malloc failed.\n",program);
      return(0);
   }
   if((locs = (float *)malloc(numpts*sizeof(float)))
      ==NULL){
      fprintf(stderr,"%s:  No dice...malloc failed.\n",program);
      return(0);
   }
   
   m1 = (double *)(image->imagedata);
   for(j=1;j<cols1+1;j++){
       for(i=0;i<numpts;i++){
           m1[(j-1)*numpts+i] = mat[j][i];
       }
   }

   for(i=0;i<numpts;i++)locs[i] = (float)mat[0][i];

   image->num_data_bands = image->row_size;
   image->row_size = numpts;
   image->col_size = 1;
   image->location = locs;
   image->location_dim = 1;
   image->location_type = VFF_LOC_EXPLICIT;

/*
** free up our junk and get the heck outa here.
*/

   free((char *)mat);
   return(1);
}

/***********************************************************************
*
*  Routine Name: 
*
*          Date:
*        
*       Purpose:  
*
*         Input: 
*
*        Output: 
*
*    Written By:  
*
* Modifications:
*
***********************************************************************/

int llde(solution,t0,tf,h,matrix,bvect,initvals,rows,cols,typea,typeb,typec)
  double **solution, *bvect, *matrix, *initvals,t0,tf,h;
  int typea,typeb,typec,rows,cols;
{
  int   idx,i,n,iscomplex=0;
  double t,h2,h6;
  complex *y,*dy=NULL,*dym=NULL, *dyt=NULL, *yt=NULL;
  unsigned numpts;
  int eval_system();

  char *program = "lmlde";

/*
** memory allocation for solution;
*/

  if(typea || typeb)iscomplex = 1;
  numpts = (int)((tf-t0)/h) + 2;
  
  for(i=0;i<cols+1;i++){
      if((solution[i] = (double *)malloc((iscomplex+1)*numpts*sizeof(double)))
         ==NULL){;
         fprintf(stderr,"%s:  [0]memory allocation failure.\n",program);
         return(0);
      }
  }

/*
** allocate memory for junk
*/

  y = (complex *)malloc((unsigned)cols * sizeof(complex));
  if(y==NULL){
     fprintf(stderr,"%s:  [1]memory allocation failure.\n",program);
     CLEANUP;
     return(0);
  }

  yt = (complex *)malloc((unsigned)cols * sizeof(complex));
  if(yt==NULL){
     fprintf(stderr,"%s:  [2]memory allocation failure.\n",program);
     CLEANUP;
     return(0);
  }

  dy = (complex *)malloc((unsigned)cols * sizeof(complex));
  if(dy==NULL){
     fprintf(stderr,"%s:  [3]memory allocation failure.\n",program);
     CLEANUP;
     return(0);
  }

  dyt = (complex *)malloc((unsigned)cols * sizeof(complex));
  if(dyt==NULL){
     fprintf(stderr,"%s:  [4]memory allocation failure.\n",program);
     CLEANUP;
     return(0);
  }

  dym = (complex *)malloc((unsigned)cols * sizeof(complex));
  if(dym==NULL){
     fprintf(stderr,"%s:  [5]memory allocation failure.\n",program);
     CLEANUP;
     return(0);
  }

/*
** set up initial values
*/

  solution[0][0] = 0.0;

  
  if(initvals!=NULL){
     for(i = 0 ; i<cols;i++){
         if(typec){
            y[i].r = solution[i+1][0] = initvals[2*i];
            y[i].i = solution[i+1][1] = initvals[2*i+1];
         }else{
            y[i].r = solution[i+1][0] = initvals[i];
            y[i].i = 0.0;
         } 
     }
  }else{
     for(i=0;i<cols;i++){
         y[i].r = y[i].i = 0.0;
     }
  }

/* 
** calculate some parameters
*/

  h2 = h/2.0;
  h6 = h/6.0;
  n  = 1;

/*
** perform runge-kutta method on the system
*/

  for(t=t0;t<tf;t+=h){
      eval_system(dy,y,t,matrix,bvect,rows,typea,typeb);

      for(i=0;i<cols;i++){
          yt[i].r = y[i].r + h2*dy[i].r;
          yt[i].i = y[i].i + h2*dy[i].i;
      }

      eval_system(dyt,yt,t+h2,matrix,bvect,rows,typea,typeb);

      for(i=0;i<cols;i++){
          yt[i].r = y[i].r + h2*dyt[i].r;
          yt[i].i = y[i].i + h2*dyt[i].i;
      }

      eval_system(dym,yt,t+h2,matrix,bvect,rows,typea,typeb);
   
      for(i=0;i<cols;i++){
          yt[i].r = y[i].r + h*dym[i].r;
          yt[i].i = y[i].i + h*dym[i].i;
          dym[i].r += dyt[i].r;
          dym[i].i += dyt[i].i;
      }

      eval_system(dyt,yt,t+h,matrix,bvect,rows,typea,typeb);
    
      for(i=0;i<cols;i++){
          y[i].r += h6*(dy[i].r + dyt[i].r + 2.0*dym[i].r);
          y[i].i += h6*(dy[i].i + dyt[i].i + 2.0*dym[i].i);
      }

      if(iscomplex){    
         idx = n*2;
         solution[0][idx] = t+h;
         for(i=0;i<cols;i++){
             solution[i+1][idx] = y[i].r; 
             solution[i+1][idx+1] = y[i].i; 
         }
      }else{
         solution[0][n] = t+h;
         for(i=0;i<cols;i++) solution[i+1][n] = y[i].r;
      }
      n++;
  }

  CLEANUP;
  return(1);
}


/***********************************************************************
*
*  Routine Name: eval_system()
*
*          Date: Sun Nov 11 20:22:38 MST 1990
*        
*       Purpose: evaluates a system represented by a matrix.  
*
*         Input: 
*
*        Output: 
*
*    Written By:  Jeremy Worley
*
* Modifications:
*
***********************************************************************/

int eval_system(dy,y,t,matrix,bvect,dim,typea,typeb)
    complex *dy, *y;
    double *matrix,*bvect,t;
    int dim,typea,typeb;
{
  int j,i;

  for(i=0;i<dim;i++){
      if(bvect!=NULL){
         if(typeb){                     
            dy[i].r = bvect[2*i];
            dy[i].i = bvect[2*i+1];
         }else{
            dy[i].r = bvect[i];
            dy[i].i = 0.0;
         }
      }else{
        dy[i].r = dy[i].i = 0.0;
      } 
      for(j=0;j<dim;j++){
          if(typea){
             dy[i].r += matrix[2*(i*dim+j)]*y[j].r-matrix[2*(i*dim+j)+1]*y[j].i;
             dy[i].i += matrix[2*(i*dim+j)]*y[j].i+matrix[2*(i*dim+j)+1]*y[j].r;
          }else{
             dy[i].r += matrix[(i*dim+j)]*y[j].r;
             dy[i].i += matrix[(i*dim+j)]*y[j].i;
          }
      }
  }

  return(1);
}
/* -library_code_end */
