/*****************************************************************************/
/* bem0_2d.f -- translated by f2c (version 19940220). */
/* File:	bem0_2d.c
 *
 * Description: constant Boundary Element solution of 2D potential problems
 *
 * Edited by:	Benny Yih
 * Modified:	03 Mar 1994
 */
/*****************************************************************************/
#ifndef lint
static const char  *sccs_id = "%W% %G% %U%%Z%";
static const char  *rcs_id = "$Header: $";
#endif
/*---------------------------------------------*/
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <limits.h>
#include <string.h>

typedef	double	real;
typedef int	integer;

/* Global variable declarations */
integer		n_boundary;
integer		n_internal;

#ifndef MAX_N_PTS
# define	MAX_N_PTS	40
#endif

/* Subroutines */
extern void	fmat_( real *, real *, real *, real *, real *, 
		       real *, real *, real *, integer *, const integer );
extern real	slnpd_( real *, real *, const integer, const integer );
extern void	inter_( real *, real *, integer *, real *, real *, 
			real *, real *, real * );
extern void	inte_( const real, const real, const real, const real,
		       const real, const real,
		       real *, real *);
extern void	inlo_( const real, const real, const real, const real,
		       real * );
extern void	input_( real *, real *, real *, real *, 
			integer *, real *,  FILE *, FILE * );
extern void	outpt_( real *, real *, real *, real *, real *,
			real *, real *, FILE * );

/*************************************************/
/*
 *     bem.f (program 1)
 *
 *     program for solution of two dimensional potential problems by the
 *     b.i.e. method with constant elements
 *
 * version:: 28 Feb 1994
 * :: 
 *
 * in:: # boundary, # internal, internal pts, boundary pts, B.C. type & value
 * out:: boundary values, interior values
 */
int main( int	argc, char *	argv[] )
{
  char *	prog_nam;

  char *	in_file = NULL;
  char *	out_file = NULL;
  FILE *	fp_in = stdin;
  FILE *	fp_out = stdout;
  extern char *	optarg;				/* getopt() flag argument */
  extern int	optind, opterr, optopt;		/* getopt() status stuff */
  int		opt_status;			/* getopt() return status */
/*---------------------------------------------*/
  /* Local variables */
  real		d;
  real		g[ MAX_N_PTS * MAX_N_PTS ];  /* was [MAX_N_PTS][MAX_N_PTS] */
  real		h[ MAX_N_PTS * MAX_N_PTS ];  /* was [MAX_N_PTS][MAX_N_PTS] */
  real		x[ MAX_N_PTS + 1 ], y[ MAX_N_PTS + 1 ];	/* boundary pt x,y */
  integer	kode[ MAX_N_PTS ];			/* B.C. type */
  real		fi[ MAX_N_PTS ];		/* B.C. value */
  real		dfi[ MAX_N_PTS ];		/* B.C. derivative value */
  real		cx[ MAX_N_PTS ], cy[ MAX_N_PTS ];	/* internal pt x,y */
  real		xm[ MAX_N_PTS ], ym[ MAX_N_PTS ];	/* boundary mid pts */
  real		sol[ MAX_N_PTS ];
  integer	nx;
/*---------------------------------------------*/
  /*     nx = maximum dimension of the system of equations. */
  nx = MAX_N_PTS;
/*---------------------------------------------*/
  prog_nam = argv[ 0 ];
  while ( ( opt_status = getopt( argc, argv, "i:o:" ) ) != EOF )
  {
    int		err_flag = 0;

    switch ( opt_status )
    {
    case 'i':					/* input file */
      in_file = optarg;
      fp_in = fopen( in_file, "r" );
      break;
    case 'o':					/* output file */
      out_file = optarg;
      fp_out = fopen( out_file, "a" );
      break;
    case '?':
      ++err_flag;
      break;
    }		/*{ opt_status }*/
    if ( err_flag )
    {
    usage:
      fprintf( stderr, "usage: %s [-i in_file] [-o out_file]\n",
	       prog_nam );
      exit( 1 );
    }		/*{ err_flag }*/
  }	/*{ opt_status != EOF }*/
/*---------------------------------------------*/

  /*     initialization of program parameters */
  input_( cx, cy, x, y, kode, fi, fp_in, fp_out );	/*     input */
  if ( fp_in != stdin )			fclose( fp_in );

  /* form system of equations */
  fmat_( x, y, xm, ym, g, h, fi, dfi, kode, nx );

  /*     solution of the system of equations */
  d = slnpd_( g, dfi, n_boundary, nx );

  /*     compute the potential values in internal points */
  inter_( fi, dfi, kode, cx, cy, x, y, sol );

  /*     output */
  outpt_( xm, ym, fi, dfi, cx, cy, sol, fp_out );

  if ( fp_out != stdout )		fclose( fp_out );
  exit( 0 );
}	/*{ main() }*/


/* n_boundary = number of coundary elements */
/* n_interal = number of internal points where the function is calculated */

void	input_( real * cx, real * cy, real * x, real * y, integer * kode,
		real * fi, FILE *	fp_in, FILE *	fp_out )
{
  integer	i;					/* Local variables */
  static char	buf[ _POSIX_PIPE_BUF ] = { '\0' };

  /* Parameter adjustments */
  --fi;
  --kode;
  --y;
  --x;
  --cy;
  --cx;

  fprintf( fp_out, " " );
  for ( i = 1; i < 80; ++i )
    fprintf( fp_out, "*" );

  if ( fgets( buf, _POSIX_PIPE_BUF, fp_in ) == NULL )		abort();

  if ( sscanf( buf, "%d,%d", & n_boundary, & n_internal ) != 2 )
    abort();

  fprintf( fp_out, "\n\n data\n\n  number of boundary elements =%3d\n  # of internal pts where the function is calculated =%3d\n",
	   n_boundary, n_internal );

  if ( n_boundary > MAX_N_PTS || n_internal > MAX_N_PTS )
  {
    fprintf( stderr, "Max # points = %d\n", MAX_N_PTS );
    exit( 1 );
  }

  /*     read internal points coordinates */
  for ( i = 1; i <= n_internal; ++i )
  {
    if ( fgets( buf, _POSIX_PIPE_BUF, fp_in ) == NULL )	abort();
    if ( sscanf( buf, "%lg,%lg", & cx[ i ], & cy[ i ] ) != 2 )
      abort();
  }		/*{ i lup }*/

  /*     read coordinates of extreme points of the boundary elements */
  /*     in array x and y */
  fprintf( fp_out, "\n\n  coordinates of the extreme points of the boundary elements\n\n    point          x                  y\n" );
    
  for ( i = 1; i <= n_boundary; ++i )
  {
    if ( fgets( buf, _POSIX_PIPE_BUF, fp_in ) == NULL )	abort();
    if ( sscanf( buf, "%lg,%lg", & x[ i ], & y[ i ] ) != 2 )
      abort();
    fprintf( fp_out, "     %3d     %14.7le     %14.7le\n",
	     i, x[ i ], y[ i ] );
  }		/*{ i lup }*/

  /*     read boundary conditions */
  /*     fi(i) = value of the potential in the node i if kode = 0, */
  /*     value of the potential derivative if kode = 1. */
  fprintf( fp_out, "\n\n  boundary conditions\n\n     node      code     prescribed value\n" );

  for ( i = 1; i <= n_boundary; ++i )
  {
    if ( fgets( buf, _POSIX_PIPE_BUF, fp_in ) == NULL )	abort();
    if ( sscanf( buf, "%d,%lg", & kode[ i ], & fi[ i ] ) != 2 )
      abort();
    fprintf( fp_out, "     %3d        %1d        %14.7le\n",
	     i, kode[ i ], fi[ i ] );
  }		/*{ i lup }*/

  return;
}	/* input_() */


/*     program 3 */
/*   this subroutine computes g and h matrices and form the system ax = f */

void	fmat_( real * x, real * y, real * xm, real * ym, real * g,
	       real * h, real * fi, real * dfi, integer * kode,
	       const integer	nx )
{
  /* System generated locals */
  integer		g_dim1, g_offset, h_dim1, h_offset;

  /* Local variables */
  integer		i, j;
  static real		ch;

  /* Parameter adjustments */
  --x;
  --y;
  --xm;
  --ym;
  --fi;
  --dfi;
  --kode;
  h_dim1 = nx;
  h_offset = h_dim1 + 1;
  h -= h_offset;
  g_dim1 = nx;
  g_offset = g_dim1 + 1;
  g -= g_offset;
#define g_ref(a_1,a_2)	g[ (a_2)*g_dim1 + (a_1) ]
#define h_ref(a_1,a_2)	h[ (a_2)*h_dim1 + (a_1) ]

  /*   compute the mid-point coordinates and store in array xm and ym */
  x[ n_boundary + 1 ] = x[ 1 ];
  y[ n_boundary + 1 ] = y[ 1 ];
  for ( i = 1; i <= n_boundary; ++i )
  {
    xm[ i ] = (x[ i ] + x[ i + 1 ]) * 0.50;
    ym[ i ] = (y[ i ] + y[ i + 1 ]) * 0.50;
  }		/*{ i lup }*/

  /*     compute g and h matrices */
  for ( i = 1; i <= n_boundary; ++i )
  {
    for ( j = 1; j <= n_boundary; ++j )
    {
      if ( i != j )
      {
	inte_( xm[ i ], ym[ i ], x[ j ], y[ j ], x[ j + 1 ], y[ j + 1 ],
	       & h_ref( i, j ), & g_ref( i, j ) );
      } else {
	inlo_( x[ j ], y[ j ], x[ j + 1 ], y[ j + 1 ], & g_ref( i, j ) );
	h_ref( i, j ) = M_PI;			/* 3.1415926 */;
      }		/*{ i != j }*/
    }		/*{ j lup }*/
  }		/*{ i lup }*/

  /*     arrange the system of equations to be ready to solve */
  for ( j = 1; j <= n_boundary; ++j )
  {
    if ( kode[ j ] > 0 )		/* derivative of potential */
      for ( i = 1; i <= n_boundary; ++i )
      {
	ch = g_ref( i, j );
	g_ref( i, j ) = - h_ref( i, j );
	h_ref( i, j ) = - ch;
      }		/*{ i lup }*/
  }		/*{ j lup }*/

  /*     dfi originally contains the independent coefficients */
  /*     after solution it will contain the values of the system unknowns */
  for ( i = 1; i <= n_boundary; ++i )
  {
    dfi[ i ] = 0.0;
    for ( j = 1; j <= n_boundary; ++j )
    {
      dfi[ i ] += h_ref( i, j ) * fi[ j ];
    }		/*{ j lup }*/
  }		/*{ i lup }*/

#undef h_ref
#undef g_ref
  return;
}		/* fmat_() */


/*     program 4 */
/*     this subroutine computes the values of the h and g matrix */
/*     off diagonal elements by means of numerical integration */
/*     along the boundary elements */

/*     dist = distance from the point under consideration to the */
/*		boundary elements */
/*     ra = distance from the point under consideration to the */
/*		integration points in the boundary elements */

void	inte_( const real	xp, const real	yp,
	       const real	x1, const real	y1,
	       const real	x2, const real	y2,
	       real *	h, real *	g )
{
  /* Local variables */
  integer	i;
  real		dist;
  real		gi[ 4 ], ax, bx, ay, by;
  real		ome[ 4 ], sig;
  real		xco[ 4 ], yco[ 4 ];

  gi[ 0 ] = 0.86113631;
  gi[ 1 ] = - gi[ 0 ];
  gi[ 2 ] = 0.33998104;
  gi[ 3 ] = - gi[ 2 ];
  ome[ 0 ] = 0.34785485;
  ome[ 1 ] = ome[ 0 ];
  ome[ 2 ] = 0.65214515;
  ome[ 3 ] = ome[ 2 ];
  ax = (x2 - x1) * 0.50;
  bx = (x2 + x1) * 0.50;
  ay = (y2 - y1) * 0.50;
  by = (y2 + y1) * 0.50;
  if ( ax != 0.0 )
  {
    register real	ta;

    ta = ay / ax;
    dist = fabs( (ta * xp - yp + y1 - ta * x1) / sqrt( ta * ta + 1) );
  } else {
    dist = fabs( xp - x1 );
  }
  sig = (x1 - xp) * (y2 - yp) - (x2 - xp) * (y1 - yp);
  if ( sig < 0.0 )
    dist = - dist;

  *g = 0.0;
  *h = 0.0;
  for ( i = 0; i < 4; ++i )
  {
    real		ra;

    xco[ i ] = ax * gi[ i ] + bx;
    yco[ i ] = ay * gi[ i ] + by;
    {
      register real	r__1, r__2;	/* System generated locals */

      r__1 = xp - xco[ i ];		/* Computing 2nd power */
      r__2 = yp - yco[ i ];		/* Computing 2nd power */
      ra = sqrt( r__1 * r__1 + r__2 * r__2 );
    }
    {
      register real	r__3;		/* System generated locals */

      r__3 = sqrt( ax * ax + ay * ay );
      *g += log( 1.0 / ra ) * ome[ i ] * r__3;
      *h -= dist * ome[ i ] * r__3 / (ra * ra);
    }
  }		/*{ i lup }*/
  return;
}		/* inte_() */

/*     program 5 */
/*     this subroutine computes the values of the diagonal */
/*     elements of the g matrix */

void	inlo_( const real	x1, const real	y1,
	       const real	x2, const real	y2,
	       real *	g )
{
  register real		ax, ay, sr;		/* Local variables */

  ax = (x2 - x1) * 0.50;
  ay = (y2 - y1) * 0.50;
  sr = sqrt( ax * ax + ay * ay );
  *g = sr * 2.0 * ( log( 1.0 / sr ) + 1.0 );
  return;
}		/* inlo_() */


/*     program 6 */
/*     solution of linear systems of equations */
/*     by the Gauss elimination method providing */
/*     for interchanging rows when encountering a zero diagonal coefficient */

/*     a : system matrix */
/*     b : originally it contains the independent coefficients. */
/*         after solution it contains the values of the system unknowns. */

/*     n0 : actual number of unknowns */
/*     nx : row and column dimension of a */

real	slnpd_( real * a, real * b, const integer n0,
		const integer nx )
{
  integer	a_dim1, a_offset;		/* System generated locals */

  /* Local variables */
  real		det;
  real		c;
  integer	i, j, k, l, k1, n1;

  /* Parameter adjustments */
  --b;
  a_dim1 = nx;
  a_offset = a_dim1 + 1;
  a -= a_offset;
#define a_ref(a_1,a_2)	a[ (a_2)*a_dim1 + (a_1) ]

  n1 = n0 - 1;
  for ( k = 1; k <= n1; ++k )
  {
    k1 = k + 1;
    c = a_ref( k, k );
    if ( fabs( c ) <= 1.0e-6 )
    {
      real	elt_abs;

      /* try to interchange rows to get non zero diagonal coefficient */
      for ( j = k1, elt_abs = 1.0e-6; j <= n0; ++j )
      {
	if (( elt_abs = fabs( a_ref( j, k ) )) > 1.0e-6 )
	{
	  for ( l = k; l <= n0; ++l )
	  {
	    c = a_ref( k, l );
	    a_ref( k, l ) = a_ref( j, l );
	    a_ref( j, l ) = c;
	  }		/*{ l lup }*/
	  c = b[ k ];
	  b[ k ] = b[ j ];
	  b[ j ] = c;
	  c = a_ref( k, k );
	  break;
	}
      }		/*{ j lup }*/
      if ( elt_abs <= 1.0e-6 )
      {
	fprintf( stderr, "*** singularity in row %4d\n", k );
	return( 0.0 );
      }
    }		/*{ fabs( c ) <= 1.0e-6 ? }*/

    /*     divide row by diagonal coefficient */
    c = a_ref( k, k );
    for ( j = k1; j <= n0; ++j )
    {
      a_ref( k, j ) /= c;
    }		/*{ j lup }*/
    b[ k ] /= c;

    /*     eliminate unknown x(k) from row i */
    for ( i = k1; i <= n0; ++i )
    {
      c = a_ref( i, k );
      for ( j = k1; j <= n0; ++j )
      {
	a_ref( i, j ) -= c * a_ref( k, j );
      }			/*{ j lup }*/
      b[ i ] -= c * b[ k ];
    }		/*{ i lup }*/
  }		/*{ k lup }*/

  /*     compute last unknown */
  if ( fabs( a_ref( n0, n0 ) ) < 1.0e-6 )
  {
    fprintf( stderr, "*** singularity in row %4d\n", k );
  }
  b[ n0 ] /= a_ref( n0, n0 );

  /*     apply backsubstitution process to compute remaining unknowns */
  for ( l = 1; l <= n1; ++l )
  {
    k = n0 - l;
    k1 = k + 1;
    for ( j = k1; j <= n0; ++j )
    {
      b[ k ] -= a_ref( k, j ) * b[ j ];
    }		/*{ j lup }*/
  }		/*{ l lup }*/

  /*     compute value of determinant */
  det = 1.0;
  for ( i = 1; i <= n0; ++i )
  {
    det *= a_ref( i, i );
  }		/*{ i lup }*/

#undef a_ref
  return( det );
}		/* slnpd_() */


/*     program 7 */
/*     this subroutine computes the potential value for internal points. */

void	inter_( real * fi, real * dfi, integer * kode, real * cx, 
		real * cy, real * x, real * y, real * sol )
{
  integer		i, j, k;		/* Local variables */

  /* Parameter adjustments */
  --sol;
  --y;
  --x;
  --cy;
  --cx;
  --kode;
  --dfi;
  --fi;

  /*     reorder fi and dfi array to put all the values of the potential */
  /*     in fi and all the values of the derivative in dfi */
  for ( i = 1; i <= n_boundary; ++i )
    if ( kode[ i ] > 0 )		/* derivative of potential */
    {
      register real	ch;

      ch = fi[ i ];
      fi[ i ] = dfi[ i ];
      dfi[ i ] = ch;
    }		/*{ kode[ i ] > 0 ? }*/

  /*     compute the potential values for internal points */
  for ( k = 1; k <= n_internal; ++k )
  {
    sol[ k ] = 0.0;
    for ( j = 1; j <= n_boundary; ++j )
    {
      real		a, b;

      inte_( cx[ k ], cy[ k ], x[ j ], y[ j ], x[ j + 1] , y[ j + 1 ],
	     & a, & b );
      sol[ k ] += dfi[ j ] * b - fi[ j ] * a;
    }		/*{ j lup }*/
    sol[ k ] /= 2.0 * M_PI;			/* 2 * 3.1415926 */
  }		/*{ k lup }*/

  return;
}		/* inter_() */

/*     program 8 */
void	outpt_( real * xm, real * ym, real * fi, real * dfi, real * cx,
		real * cy, real * sol, FILE * fp_out )
{
  integer	i, k;			/* Local variables */

  /* Parameter adjustments */
  --sol;
  --cy;
  --cx;
  --dfi;
  --fi;
  --ym;
  --xm;

  fprintf( fp_out, " " );
  for ( i = 1; i < 80; ++i )
    fprintf( fp_out, "*" );
  fprintf( fp_out, "\n\n results\n\n  boundary nodes\n\n         x                y            potential          potential derivative\n" );

  for ( i = 1; i <= n_boundary; ++i )
  {
    fprintf( fp_out, "  %14.7le      %14.7e      %14.7le      %14.7le\n",
	     xm[ i ], ym[ i ], fi[ i ], dfi[ i ] );
  }		/*{ i lup }*/

  fprintf( fp_out, "\n\n  internal points\n\n           x                  y              potential\n\n" );
  for ( k = 1; k <= n_internal; ++k )
  {
    fprintf( fp_out, "     %14.7le     %14.7e     %14.7le\n",
	     cx[ k ], cy[ k ], sol[ k ] );
  }

  fprintf( fp_out, " " );
  for ( i = 1; i < 80; ++i )
    fprintf( fp_out, "*" );
  fprintf( fp_out, "\n" );

  return;
}		/* outpt_() */