/********************************* hs.c ************************************
                          Xinyu Deng (cindy)
                             Feb.26, 1994
***************************************************************************/
/************************************************************************
have made modification: eliminate extra 2 multiple for B
allow dynamic time step
add output relative arc length and velocity to output2
*************************************************************************/

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

/* constants definition */
#define MAX_LEN  80
#define SEPCHARS " !=#\n"

#ifndef PI
#define PI       3.141592653589793238462643383280
#endif

#ifndef ERROR
#define ERROR    -1   
#endif

#ifndef TRUE
#define TRUE    1   
#endif

#ifndef FALSE
#define FALSE    0   
#endif

/* function prototype */
extern void allocate_memory(void);
extern void free_memory(void);
extern double **allocate_2D_double_array(int row,int col);
extern void free_2D_double_array(double **pointer);
extern void readinit(FILE *input);
extern int cm_Ri(int);
extern int cm_Li(int);
extern void cm_kg(double*, double*);
extern void cm_p(void);
extern int hs_iteration(void);

extern int lu_decomp(double** a,int n,int* indx,int* d);
extern void lu_bksb(double** a, int n, int* indx, double* b);
extern int lu_lin_solver(double** a, int n, double* b, double* x);
extern int lu_invmat(double** mat, int n, double** invmat);
extern int lu_det(double** mat, int n, double* det);

/*** global variable declaration ***/
double delta_t, delta1, time;
double beta;
int k, n, m, mn;
double area, s;
FILE *output1, *output2, *output3, *output4;

double *x, *y;        /* x, y coordinate of each point */
double *zx, *zy;        /* x, y coordinate of each point */
double *tRx, *tRy;
double *tLx, *tLy;
double *nLx, *nLy;
double *nRx, *nRy;
double *tRLx, *tRLy;
double *dR, *dL, *dRL, *d;
double *tox, *toy;
double *normalx, *normaly;
double *kv;           /* curvature */
double *u;            /* mn+1 unknowns: v1...vn, c */
double *b, *bl, *br;
double *p;
double *f;
double *rel_s;
double *abs_s;
double **a, **g, **ginv;           /* mn+1 by mn+1 matrix */

/*---------------------------------------------------------------------------
        main interface
----------------------------------------------------------------------------*/
void main(int argc, char *argv[])
{
  int i;
  FILE *input;
  char *prog = argv[0];
  char buf[100];
  char *str;

/* check if both input file and output file are given */
  if( argc < 6 )
  {
        fprintf(stderr, "Usage: %s <in~> <re~> <gre~> <n~> <f~>\n", prog);
        exit( 1 );
  }

/* open input file */
  if( (input = fopen(argv[1], "r")) == NULL )
  {
        fprintf(stderr, "%s: can't open %s\n", prog, argv[1]);
        exit( 1 );
  }

/* read parameters from input file */
  readinit(input);
  fclose( input );
/* input k and delta_t */
  printf("\nEnter number of iterations k: ");
  fflush( stdin );
  scanf(" %d", &k);
  if( k < 1 )
  {
        fprintf(stderr, "Invalid number of iteration: %d\n", k);
        exit( -1 );
  }
  printf("\nEnter time step delta_t: ");
  fflush( stdin);
  gets( buf );
  str = strtok( buf, SEPCHARS);
  delta_t = atof(str);
  printf("\nEnter beta: ");
  fflush( stdin);
  gets( buf );
  str = strtok( buf, SEPCHARS);
  beta = atof(str);
  printf("\n");
  printf("m=%d n=%d k=%d beta=%lf  delta_t=%lf\n", m, n, k, beta, delta_t);

/* open output file1 */
  if( (output1 = fopen(argv[2], "w")) == NULL )
  {
    printf("%s: can't open %s for output\n", prog, argv[2]);
    exit( 1 );
  }
  fprintf(output1, "m=%d n=%d   k=%d  beta=%lf  delta_t=%lf\n", m, n, k, beta, delta_t);

/* open output file2 */
  if( (output2 = fopen(argv[3], "w")) == NULL )
  {
    printf("%s: can't open %s for output\n", prog, argv[3]);
    exit( 1 );
  }
  fprintf(output2, "#m=%d n=%d   k=%d  beta=%lf  delta_t=%lf\n", m, n, k, beta, delta_t);

/* open output file3 */
  if( (output3 = fopen(argv[4], "w")) == NULL )
  {
    printf("%s: can't open %s for output\n", prog, argv[4]);
    exit( 1 );
  }
  fprintf(output3, "#k=%d   m=%d   n=%d   beta=%lf  delta_t=%lf\n", k, m, n, beta, delta_t);

/* open output file4 */
  if( (output4 = fopen(argv[5], "w")) == NULL )
  {
    printf("%s: can't open %s for output\n", prog, argv[5]);
    exit( 1 );
  }
  fprintf(output4, "#m=%d n=%d   k=%d  beta=%lf  delta_t=%lf\n", m, n, k, beta, delta_t);


/* do calculation here */
  if( hs_iteration() == ERROR ) {
    printf("singular matrix found in hs_iteration\n");
  }
  fclose( output1 );
  fclose( output2 );
  fclose( output3 );
  fclose( output4 );
/* free memory */
  free_memory();
}  /* end of main */

/*---------------------------------------------------------------------------
         cm_Ri() - this function calculate Ri
--------------------------------------------------------------------------*/
int cm_Ri(int i)
{
  int Ri;

  if (i == i/n*n)       /* i/n = integer */
    Ri = i+1-n;
  else 
    Ri = i+1;
  return Ri;
}

/*---------------------------------------------------------------------------
         cm_Li() - this function calculate Ri
--------------------------------------------------------------------------*/
int cm_Li(int i)
{
  int Li;


  if (i-1 == (i-1)/n*n ) /* (i-1)/n = integer */
    Li = i-1+n;
  else 
    Li = i-1;
  return Li;
}

/*---------------------------------------------------------------------------
         readinit() - this function read initial value from input file
---------------------------------------------------------------------------*/
void readinit(FILE *input)
{
  char line[MAX_LEN];
  char *tok;
  int i=1;

  while( fgets(line, MAX_LEN, input) ) {
    if (line[0] == 'm' ) {
      tok = strchr(line, 'm');
      tok++;
      tok = strtok(tok, SEPCHARS);
      m = atoi( tok );
    } else if (line[0] == 'n' ) {
      tok = strchr(line, 'n');
      tok++;
      tok = strtok(tok, SEPCHARS);
      n = atoi( tok );
      mn = m*n;
/* allocate memory for all variables */
      allocate_memory();
    } else if( line[0] == ' ') {
      tok = strtok(line, SEPCHARS);
      if( ! tok ) continue;
      x[i] = atof( tok );
      tok = strtok(NULL, SEPCHARS);
      y[i] = atof( tok );
      i++;
    }
  } /* end of while */
} /* end of readinit() */



/*---------------------------------------------------------------------------
         cm_kg() - this function calculate g matrix
---------------------------------------------------------------------------*/
void cm_kg(double *x1, double *y1)

{
  int i, j, Ri, Li, Rj, Lj;
  double dp, dRj, dLj, dj, dij;
  double a1, a2, a3, a4 ;
  double tjix,tjiy,A,B;

  for (i=1; i<=mn; i++)
  {
    Ri = cm_Ri(i);
    Li = cm_Li(i);

    tRx[i] = x1[Ri] - x1[i];
    tRy[i] = y1[Ri] - y1[i];
    tLx[i] = x1[i] - x1[Li];
    tLy[i] = y1[i] - y1[Li];
    tRLx[i] = x1[Ri] - x1[Li];
    tRLy[i] = y1[Ri] - y1[Li];

    dR[i] = sqrt(tRx[i]*tRx[i]+tRy[i]*tRy[i]);
    dL[i] = sqrt(tLx[i]*tLx[i]+tLy[i]*tLy[i]);
    dRL[i] = sqrt(tRLx[i]*tRLx[i]+tRLy[i]*tRLy[i]);

    tLx[i]=tLx[i]/dL[i];
    tRx[i]=tRx[i]/dR[i];
    tLy[i]=tLy[i]/dL[i];
    tRy[i]=tRy[i]/dR[i];
      
    nRx[i] = tRy[i];
    nRy[i] = -tRx[i];
    nLx[i] = tLy[i];
    nLy[i] = -tLx[i];

    d[i] = 0.5*(dR[i]+dL[i]);
    tox[i] = (dR[i]*tLx[i]+dL[i]*tRx[i])/dRL[i];
    toy[i] = (dR[i]*tLy[i]+dL[i]*tRy[i])/dRL[i];
    normalx[i] = toy[i];
    normaly[i] = -tox[i];
    kv[i] = -2.0*(tRx[i]*nLx[i]+tRy[i]*nLy[i])/dRL[i];
/*    printf("i= %d  kv= %lf\n", i,kv[i]); */
  }
   
  /* computes vector bl, b, br */
  for (i=1; i<=mn; i++)
  {
    Ri = cm_Ri(i);
    Li = cm_Li(i);
    bl[i] =  delta1*(-2*(normalx[Li]*nLx[i]+normaly[Li]*nLy[i])*
                     (nRx[i]*nLx[i]+nRy[i]*nLy[i])/(dL[i]*dRL[i])
                    +
                     kv[i]*(tRLx[i]*normalx[Li]+tRLy[i]*normaly[Li])
                     /(dRL[i]*dRL[i]));
    b[i] =   delta1*( 2*(normalx[i]*nLx[i]+normaly[i]*nLy[i])*
                     (nRx[i]*nLx[i]+nRy[i]*nLy[i])/(dL[i]*dRL[i])   
                    +
                     2*(normalx[i]*nRx[i]+normaly[i]*nRy[i])*
                     (nRx[i]*nLx[i]+nRy[i]*nLy[i])/(dR[i]*dRL[i]));
    br[i] =  delta1*( -2*(normalx[Ri]*nRx[i]+normaly[Ri]*nRy[i])*
                     (nRx[i]*nLx[i]+nRy[i]*nLy[i])/(dR[i]*dRL[i])
                    -
                     kv[i]*(tRLx[i]*normalx[Ri]+tRLy[i]*normaly[Ri]) 
                     /(dRL[i]*dRL[i]));
  }

  /* computes matrix a */
  for (i=1; i<=mn; i++)
  {
    Ri = cm_Ri(i);
    Li = cm_Li(i);
    for (j=1; j<=mn; j++)
    {
      Rj = cm_Ri(j);
      Lj = cm_Li(j);
      dij = sqrt((x1[i]-x1[j])*(x1[i]-x1[j])+(y1[i]-y1[j])*(y1[i]-y1[j]));
 
      if (i==j)
            {   tjix=tox[i];
               tjiy = toy[i];}
       else
            {   tjix= (x1[j]-x1[i])/dij;
                tjiy= (y1[j]-y1[i])/dij; }

       A= -tox[i]*tjix - toy[i]*tjiy;
       B= tox[j]*tjix + toy[j]*tjiy;

      if ( A<0.01 && A > 0.0 )  A=0.01;
      else if ( A>-0.01 && A<=0.0) A=-0.01;

      if ( B<0.01 && B>=0.0)  B= 0.01;
      else if ( B>-0.01 && B<0.0) B=-0.01;

      
     
             
               
      a1 = dij + 0.5*A*dR[i] +0.5*B*dR[j]; 
      a2 = dij - 0.5*A*dL[i] -0.5*B*dL[j]; 
      a3 = dij + 0.5*A*dR[i] -0.5*B*dL[j]; 
      a4 = dij - 0.5*A*dL[i] +0.5*B*dR[j]; 

      a[j][i] = ( a1*a1*(log(fabs(a1)+0.00000001)-1.5)
                 +a2*a2*(log(fabs(a2)+0.00000001)-1.5)
                 -a3*a3*(log(fabs(a3)+0.00000001)-1.5)
                 -a4*a4*(log(fabs(a4)+0.00000001)-1.5) )
                /A/B/(dR[i]+dL[i]); 
      a[j][i] = a[j][i]/2.0/PI;
    }  /* end for j */
  }   /* end for i */
/* 
  for (j=1; j<=m; j++)
  {
    for (i=j*n+1; i<=(j+1)*n; i++)
    {
      sum_a += a[i][j*n+1];
      
    }
  }
*/

  for (i=1; i<=mn; i++)
  {
    Ri = cm_Ri(i);
    Li = cm_Li(i);
    a[Li][i] = a[Li][i] -beta*bl[i];
    a[i][i] = a[i][i] - beta*b[i];
    a[Ri][i] = a[Ri][i] -beta* br[i];
  }

  /* computes matrix g */
  g[mn][mn] = 0.0;
  for (i=0; i<mn; i++)
  {
 /*   printf("i= %d  g= %lf\n", i, g[i][i]); */
    g[mn][i] = 1.0;                    /* last column */
    g[i][mn] = d[i+1];      /* last row */
    for (j=0; j<mn; j++)

    {
      g[j][i] = a[j+1][i+1];
    }
  }
} /* end of cm_kg() */

/*---------------------------------------------------------------------------
        cm_p() - this function  calculate vector p
---------------------------------------------------------------------------*/
void cm_p()
{
  int i;

  for (i=0; i<mn; i++)
  {
    p[i] = kv[i+1];
  }
  p[mn] = 0.0;
}

/*---------------------------------------------------------------------------
        hs_iteration() - this function performs calculations for
                         heleshow problem
---------------------------------------------------------------------------*/
int hs_iteration()
{
  int cnt, i, j, i1, Ri, Li, Rj, Lj, i2, j2;
  double maxu, lamda;
  double timeprint;
  double kavg, kag, tmps;
  int Disapr = FALSE;

  printf("mn= %d\n", mn);

  time = 0.0;
  timeprint = 0.0;
 delta1= delta_t*beta;
/*
 delta1= 0.001;
 delta1= 0.039905;
 delta1= 0.032809;
*/
  for ( cnt=0; cnt<=k; cnt++)

  {
    cm_kg(x,y); 
    /* calculate area and s */
    area = 0.0;
    s = 0.0;
    for(i=1; i<=mn; i++)
    {
      Ri = cm_Ri(i);
      Li = cm_Li(i);
      area += 0.5*x[i]*(y[Ri]-y[Li]); 
/*
      area += 0.5*fabs(y[i]*(x[Ri]-x[Li]));
*/
      s += dR[i];
    }

    /* compute relative arc length rel_s */
    tmps= 0.0;
    abs_s[1] = 0.0;
    rel_s[1] = 0.0;
    for(i=2; i<=mn; i++)
    {
      tmps += dL[i];
      abs_s[i] = tmps;
      rel_s[i] = abs_s[i]/s;
    }

   /* check disappearence */
    for(j=0; j<=m-1; j++)
    {
      kavg = 0.0;
      for(i=1; i<=n; i++)
      {
        kavg += fabs(kv[j*n+i]);
      }
      kavg = kavg/(double)n;

      if (kavg > 200)
      {
        Disapr = TRUE;
        printf("one piece vanish !\n");

        /* output data to output4 file */
        fprintf(output4, "#%d  t=%lf  delta_t=%lf\n", cnt-1, time-delta1, delta1);
        fprintf(output4, "#a=%20.15f, s=%20.15f\n", area, s);
        for (j2=0; j2<m; j2++)
        {
          for (i2=j2*n+1; i2<=(j2+1)*n; i2++)
          {
            fprintf(output4, " %20.15f   %20.15f\n", x[i2], y[i2]);
          }
          fprintf(output4, "%20.15f   %20.15f\n", x[j2*n+1], y[j2*n+1]);
          fprintf(output4, "\n");
        }

        if (j<m)
        {
          for( i1=j*n+1; i1<=mn-n; i1++)
          {
            x[i1] = x[i1+n];
            y[i1] = y[i1+n];
            kv[i1] = kv[i1+n];
          }
        }
        m = m-1;
        mn = m*n;
      }
    }

    /* output time, area, s, kv in <output file3>*/
    kag = 0.0;
    for(i=1; i<=n; i++)
    {
      kag += fabs(kv[i]);
    }
    kag = kag/(double)n;
    fprintf(output3, " %20.15f   %20.15f    %20.15f   %20.15f \n", time, kag, area, s);

   /* computes vector p */
    cm_p();

   /* LU decomposition to solve for u */
    if ( lu_lin_solver(g, mn+1, p, u) == ERROR ) {
      return (ERROR);
    }
    
    /* output x,y,kv, area, s*/
    if ( time >= timeprint || cnt == k || Disapr == TRUE)
    {
      fprintf(output2, "#%d  t=%lf  delta_t=%lf\n", cnt, time, delta1);
      fprintf(output2, "#a=%20.15f, s=%20.15f\n", area, s);
      fprintf(output1, "#%d  t=%lf  delta_t=%lf\n", cnt, time, delta1);
      fprintf(output1, "a=%20.15f, s=%20.15f\n", area, s);


      for (j=0; j<m; j++)
      {
        fprintf(output1, "[\n");
        for (i=j*n+1; i<=(j+1)*n; i++)
        {
          fprintf(output2, " %20.15f   %20.15f  %20.15f  %20.15f  %20.15f %20.15f\n", x[i], y[i], abs_s[i], rel_s[i], kv[i], u[i-1]);
          fprintf(output1, " %20.15f   %20.15f   %20.15f   %20.15f\n", x[i], y[i
], kv[i], u[i-1]);
        }
        fprintf(output2, "%20.15f   %20.15f  %20.15f  %20.15f %20.15f %20.15f\n", x[j*n+1], y[j*n+1], s, 1.0, kv[j*n+1], u[j*n]);
        fprintf(output2, "\n");
        fprintf(output1, "!%20.15f   %20.15f   %20.15f   %20.15f\n",
               x[j*n+1], y[j*n+1], kv[j*n+1], u[j*n]);
        fprintf(output1, "]\n");
      } /* end of for */
      timeprint += k*delta_t*0.05;
      Disapr = FALSE;
   } /* end of if */

    maxu = 0.0; 
    for(j=0; j < mn; j++ )
    {
       if (maxu < fabs(u[j]) )
          maxu = fabs(u[j]);
    }
    /* synamic time step */
    if (maxu > 1)
      delta1 = 1.0/maxu * delta_t;  
    else delta1=delta_t;
    printf("time step= %d,  delta_t= %lf\n", cnt, delta1);

   /* update new x, y */
    for(i=1; i<=mn; i++)
    {
      x[i] = x[i] + delta1*normalx[i]*u[i-1];
      y[i] = y[i] + delta1*normaly[i]*u[i-1];
    }
    
   /* rearrangement of points */
    for (i=1; i<=mn; i++)
    {
     Ri = cm_Ri(i);
     Li = cm_Li(i);

     tRx[i] = x[Ri] - x[i];
     tRy[i] = y[Ri] - y[i];
     tLx[i] = x[i] - x[Li];
     tLy[i] = y[i] - y[Li];

     dR[i] = sqrt(tRx[i]*tRx[i]+tRy[i]*tRy[i]);
     dL[i] = sqrt(tLx[i]*tLx[i]+tLy[i]*tLy[i]);

     tLx[i]=tLx[i]/dL[i];
     tRx[i]=tRx[i]/dR[i];
     tLy[i]=tLy[i]/dL[i];
     tRy[i]=tRy[i]/dR[i];
    }

    for (i=1; i<=mn; i++)
    {
      lamda = (dR[i]-dL[i])/2.0/(1.0+tLx[i]*tRx[i]+tLy[i]*tRy[i]);
      x[i] = x[i] + lamda*(tLx[i]+tRx[i]);
      y[i] = y[i] + lamda*(tLy[i]+tRy[i]);
    }

  /* end of rearrangement */


    printf("time= %lf\n", time);
    time += delta1;
  } /* end of for */

} /* end of hs_iteration()



/*---------------------------------------------------------------------------
        allocate_memory() - this function allocates memory for all
                        variables.
----------------------------------------------------------------------------*/
void allocate_memory()
{
        x = (double *) calloc(mn+1, sizeof( double ) );
        if ( ! x )
        {
           fprintf(stderr, "Out of memory...\n");
           exit ( -1 );
        }

        y = (double *) calloc(mn+1, sizeof( double ) );
        if ( ! y )
        {
           fprintf(stderr, "Out of memory...\n");
           exit ( -1 );
        }

        zx = (double *) calloc(mn+1, sizeof( double ) );
        if ( ! zx )
        {
           fprintf(stderr, "Out of memory...\n");
           exit ( -1 );
        }

        zy = (double *) calloc(mn+1, sizeof( double ) );
        if ( ! zy )
        {
           fprintf(stderr, "Out of memory...\n");
           exit ( -1 );
        }

        tRx = (double *) calloc(mn+1, sizeof( double ) );
        if ( ! tRx )
        {
           fprintf(stderr, "Out of memory...\n");
           exit ( -1 );
        }

        tRy = (double *) calloc(mn+1, sizeof( double ) );
        if ( ! tRy )
        {
           fprintf(stderr, "Out of memory...\n");
           exit ( -1 );
        }

        tLx = (double *) calloc(mn+1, sizeof( double ) );
        if ( ! tLx )
        {
           fprintf(stderr, "Out of memory...\n");
           exit ( -1 );
        }

        tLy = (double *) calloc(mn+1, sizeof( double ) );
        if ( ! tLy )
        {
           fprintf(stderr, "Out of memory...\n");
           exit ( -1 );
        }

        nLx = (double *) calloc(mn+1, sizeof( double ) );
        if ( ! nLx )
        {
           fprintf(stderr, "Out of memory...\n");
           exit ( -1 );
        }

        nLy = (double *) calloc(mn+1, sizeof( double ) );
        if ( ! nLy )
        {
           fprintf(stderr, "Out of memory...\n");
           exit ( -1 );
        }

        nRx = (double *) calloc(mn+1, sizeof( double ) );
        if ( ! nRx )
        {
           fprintf(stderr, "Out of memory...\n");
           exit ( -1 );
        }

        nRy = (double *) calloc(mn+1, sizeof( double ) );
        if ( ! nRy )
        {
           fprintf(stderr, "Out of memory...\n");
           exit ( -1 );
        }

        tRLx = (double *) calloc(mn+1, sizeof( double ) );
        if ( ! tRLx )
        {
           fprintf(stderr, "Out of memory...\n");
           exit ( -1 );
        }

        tRLy = (double *) calloc(mn+1, sizeof( double ) );
        if ( ! tRLy )
        {
           fprintf(stderr, "Out of memory...\n");
           exit ( -1 );
        }

        dR = (double *) calloc(mn+1, sizeof( double ) );
        if ( ! dR )
        {
           fprintf(stderr, "Out of memory...\n");
           exit ( -1 );
        }

        dL = (double *) calloc(mn+1, sizeof( double ) );
        if ( ! dL )
        {
           fprintf(stderr, "Out of memory...\n");
           exit ( -1 );
        }

        dRL = (double *) calloc(mn+1, sizeof( double ) );
        if ( ! dRL )
        {
           fprintf(stderr, "Out of memory...\n");
           exit ( -1 );
        }

        d = (double *) calloc(mn+1, sizeof( double ) );
        if ( ! d )
        {
           fprintf(stderr, "Out of memory...\n");
           exit ( -1 );
        }

        tox = (double *) calloc(mn+1, sizeof( double ) );
        if ( ! tox )
        {
           fprintf(stderr, "Out of memory...\n");
           exit ( -1 );
        }

        toy = (double *) calloc(mn+1, sizeof( double ) );
        if ( ! toy )
        {
           fprintf(stderr, "Out of memory...\n");
           exit ( -1 );
        }

        normalx = (double *) calloc(mn+1, sizeof( double ) );
        if ( ! normalx )
        {
           fprintf(stderr, "Out of memory...\n");
           exit ( -1 );
        }

        normaly = (double *) calloc(mn+1, sizeof( double ) );
        if ( ! normaly )
        {
           fprintf(stderr, "Out of memory...\n");
           exit ( -1 );
        }

        kv = (double *) calloc(mn+1, sizeof( double ) );
        if ( ! kv )
        {
           fprintf(stderr, "Out of memory...\n");
           exit ( -1 );
        }

        u = (double *) calloc(mn+1, sizeof( double ) );
        if ( ! u )
        {
           fprintf(stderr, "Out of memory...\n");
           exit ( -1 );
        }

        bl = (double *) calloc(mn+1, sizeof( double ) );
        if ( ! bl )
        {
           fprintf(stderr, "Out of memory...\n");
           exit ( -1 );
        }

        b = (double *) calloc(mn+1, sizeof( double ) );
        if ( ! b )
        {
           fprintf(stderr, "Out of memory...\n");
           exit ( -1 );
        }

        br = (double *) calloc(mn+1, sizeof( double ) );
        if ( ! br )
        {
           fprintf(stderr, "Out of memory...\n");
           exit ( -1 );
        }

        p = (double *) calloc(mn+1, sizeof( double ) );
        if ( ! p )
        {
           fprintf(stderr, "Out of memory...\n");
           exit ( -1 );
        }

        f = (double *) calloc(mn+1, sizeof( double ) );
        if ( ! f )
        {
           fprintf(stderr, "Out of memory...\n");
           exit ( -1 );
        }

        abs_s = (double *) calloc(mn+1, sizeof( double ) );
        if ( ! abs_s )
        {
           fprintf(stderr, "Out of memory...\n");
           exit ( -1 );
        }

        rel_s = (double *) calloc(mn+1, sizeof( double ) );
        if ( ! rel_s )
        {
           fprintf(stderr, "Out of memory...\n");
           exit ( -1 );
        }

        if ( (a = allocate_2D_double_array(mn+1, mn+1)) == NULL )
        {
           fprintf(stderr, "Out of memory...\n");
           exit ( -1 );
        }

        if ( (g = allocate_2D_double_array(mn+1, mn+1)) == NULL )
        {
           fprintf(stderr, "Out of memory...\n");
           exit ( -1 );
        }

        if ( (ginv = allocate_2D_double_array(mn+1, mn+1)) == NULL )
        {
           fprintf(stderr, "Out of memory...\n");
           exit ( -1 );
        }
}   /* end of allocate_memory */

/*---------------------------------------------------------------------------
        free_memory - this function is used to free memory for all
                variables.
----------------------------------------------------------------------------*/
void free_memory()
{
        if( x ) free ( ( void * ) x );
        if( y ) free ( ( void * ) y );
        if( zx ) free ( ( void * ) zx );
        if( zy ) free ( ( void * ) zy );
        if( tRx ) free ( ( void * ) tRx );
        if( tRy ) free ( ( void * ) tRy );
        if( tLx ) free ( ( void * ) tLx );
        if( tLy ) free ( ( void * ) tLy );
        if( nRx ) free ( ( void * ) nRx );
        if( nRy ) free ( ( void * ) nRy );
        if( nLx ) free ( ( void * ) nLx );
        if( nLy ) free ( ( void * ) nLy );
        if( tRLx ) free ( ( void * ) tRLx );
        if( tRLy ) free ( ( void * ) tRLy );
        if( dR ) free ( ( void * ) dR );
        if( dL ) free ( ( void * ) dL );
        if( dRL ) free ( ( void * ) dRL );
        if( d ) free ( ( void * ) d );
        if( tox ) free ( ( void * ) tox );
        if( toy ) free ( ( void * ) toy );
        if( normalx ) free ( ( void * ) normalx );
        if( normaly ) free ( ( void * ) normaly );
        if( kv ) free ( ( void * ) kv );
        if( u ) free ( ( void * ) u );
        if( bl ) free ( ( void * ) bl );
        if( br ) free ( ( void * ) br );
        if( b ) free ( ( void * ) b );
        if( p ) free ( ( void * ) p );
        if( f ) free ( ( void * ) f );
        if( rel_s ) free ( ( void * ) rel_s );
        if( abs_s ) free ( ( void * ) abs_s );
        free_2D_double_array( a );
        free_2D_double_array( g );
        free_2D_double_array( ginv );
}   /* end of free_memory */

double **allocate_2D_double_array(int row,int col)
{
        int i;
        double **prow, *pdata;

        pdata=(double *) calloc(1, row*col*sizeof( double ));
        if(pdata ==  NULL) return (NULL);
        prow=(double **) calloc(1, row*sizeof(double *));
        if(prow ==  NULL) return (NULL);

        for(i=0; i<row; i++)
        {
           prow[i] = pdata;
           pdata += col;
        }

        return (prow);
}

void free_2D_double_array(double **pointer)
{
        if ( pointer && * pointer )
        {
           free((void *) * pointer);
           free((void *) pointer);
        }
        else if( pointer )
        {
           free((void *) pointer);
        }
}