/*----------------------------------------------------------------------*/
/*									*/
/*  baff								*/
/*	Compute the location and size of baffles within a length of	*/
/*  tubing.								*/
/*									*/
/*  The input arguments are determined using getopt.  Floating point	*/
/*  conversion of dimensional arguments is performed using strtod.  No	*/
/*  other special syntax checking is performed.				*/
/*									*/
/*  The reflectance point is computed by considering the similar	*/
/*  triangles that occur on either side of it extending up to the	*/
/*  extreme point of the entry and exit apertures, across these		*/
/*  apertures to the opposite tube wall and back to the reflectance	*/
/*  point.								*/
/*									*/
/*  The baffles in front of the reflectance point are computed first.	*/
/*  These are found by finding the intersection of the line from the	*/
/*  extreme edge of the entry aperture to the leading shadow point	*/
/*  with the line joining the proximal entry and exit aperture points.	*/
/*  Allowance here is made for both baffle tolerance and shadow		*/
/*  overlap.								*/
/*									*/
/*  A new leading shadow point is then computed on the tube wall using	*/
/*  the just computed baffle position and the extreme exit aperture	*/
/*  position.  The baffle tolerance is used here to produce a minimal	*/
/*  shadow coverage position.						*/
/*									*/
/*  This same process is then repeated for the baffles following the	*/
/*  reflectance point, reversing the role of entry and exit apertures.	*/
/*									*/
/*----------------------------------------------------------------------*/

#include <stdlib.h>
#include <stdio.h>
#include <unistd.h>

static double min_tube_diameter   =  200.0;
static double tube_length         = 1000.0;
static double entry_semidiameter  =   75.0;
static double exit_semidiameter   =   75.0;
static double shadow_overlap      =    1.0;
static double baffle_tolerance    =    1.0;

#define MAX_BAFFLES 128
#define MIN_REMAINING 0.000001

static double tube_semidiameter;
static double reflection_z;
static double lead_shadow_z;
static double trail_shadow_z;

static int    num_lead_baffles;
static double lead_z[ MAX_BAFFLES ];
static double lead_y[ MAX_BAFFLES ];

static int    num_trail_baffles;
static double trail_z[ MAX_BAFFLES ];
static double trail_y[ MAX_BAFFLES ];

static double reflection_point_z(
  const double yf,
  const double y0,
  const double z0,
  const double y1,
  const double z1 );

static double intersection_z(
  const double y0a,
  const double z0a,
  const double y1a,
  const double z1a,
  const double y0b,
  const double z0b,
  const double y1b,
  const double z1b );

static double baffle_semidiameter_y(
  const double zf,
  const double y0,
  const double z0,
  const double y1,
  const double z1 );

static double lead_point_z(
  const double yf,
  const double y0,
  const double z0,
  const double y1,
  const double z1 );

main( int argc, char **argv )
    {
    int baffle_index;
    double previous_z;
    int opt_char;

    while( ( opt_char = getopt( argc, argv, "b:d:e:l:s:x:" ) ) != EOF )
	{
	switch( opt_char )
	    {
	case 'd':
	    min_tube_diameter = strtod( optarg, NULL );
	    break;
	case 'l':
	    tube_length = strtod( optarg, NULL );
	    break;
	case 'e':
	    entry_semidiameter = strtod( optarg, NULL );
	    break;
	case 'x':
	    exit_semidiameter = strtod( optarg, NULL );
	    break;
	case 's':
	    shadow_overlap = strtod( optarg, NULL );
	    break;
	case 'b':
	    baffle_tolerance = strtod( optarg, NULL );
	    break;
	    }
	}

    /* The coordinate system here is one for which the	*/
    /* z value is zero at the entry.  The y value is	*/
    /* zero along the symmetry axis. 			*/

    /* Compute the position of the reflection point.	*/

    tube_semidiameter = 0.5 * min_tube_diameter;
    reflection_z =
      reflection_point_z(
	tube_semidiameter,
	- entry_semidiameter,
	0.0,
	- exit_semidiameter,
	tube_length );

    /* Find the baffles that precede the reflection	*/
    /* point.						*/

    baffle_index = 0;
    lead_shadow_z = reflection_z - 0.5 * shadow_overlap;
    while(
      lead_shadow_z > MIN_REMAINING &&
      baffle_index < MAX_BAFFLES )
	{
	/* Find the baffle that, given the tolerance on	*/
	/* the baffles, will overlap the shadow by the	*/
	/* given amount.				*/

	lead_z[ baffle_index ] =
	  intersection_z(
	    - entry_semidiameter,
            0.0,
	    tube_semidiameter,
	    lead_shadow_z + shadow_overlap,
	    entry_semidiameter + baffle_tolerance,
	    0.0,
	    exit_semidiameter + baffle_tolerance,
	    tube_length );

	/* Find the design baffle semi-diameter at this	*/
	/* point.					*/

	lead_y[ baffle_index ] =
	    baffle_semidiameter_y(
	      lead_z[ baffle_index ],
	      entry_semidiameter,
	      0.0,
	      exit_semidiameter,
	      tube_length );
	      
	/* Find the new lead shadow point given the	*/
	/* baffle tolerance.				*/

	lead_shadow_z =
	  lead_point_z(
	    tube_semidiameter,
	    lead_y[ baffle_index ] + baffle_tolerance,
	    lead_z[ baffle_index ],
	    - exit_semidiameter,
	    tube_length );
	
	++baffle_index;
	}
    num_lead_baffles = baffle_index;

    /* Find the baffles that follow the reflection	*/
    /* point.						*/

    baffle_index = 0;
    trail_shadow_z = reflection_z + 0.5 * shadow_overlap;
    while(
      trail_shadow_z < tube_length - MIN_REMAINING &&
      baffle_index < MAX_BAFFLES )
	{
	/* Find the baffle that, given the tolerance on	*/
	/* the baffles, will overlap the shadow by the	*/
	/* given amount.				*/

	trail_z[ baffle_index ] =
	  intersection_z(
	    - exit_semidiameter,
            tube_length,
	    tube_semidiameter,
	    trail_shadow_z - shadow_overlap,
	    entry_semidiameter + baffle_tolerance,
	    0.0,
	    exit_semidiameter + baffle_tolerance,
	    tube_length );

	/* Find the design baffle semi-diameter at this	*/
	/* point.					*/

	trail_y[ baffle_index ] =
	    baffle_semidiameter_y(
	      trail_z[ baffle_index ],
	      entry_semidiameter,
	      0.0,
	      exit_semidiameter,
	      tube_length );
	      
	/* Find the new lead shadow point given the	*/
	/* baffle tolerance.				*/

	trail_shadow_z =
	  lead_point_z(
	    tube_semidiameter,
	    trail_y[ baffle_index ] + baffle_tolerance,
	    trail_z[ baffle_index ],
	    - entry_semidiameter,
	    0.0 );
	
	++baffle_index;
	}
    num_trail_baffles = baffle_index;

    /* Print the results.				*/

    (void)fprintf(
      stdout,
      "    %15s %15s %15s\n",
      "Separation",
      "Position",
      "Semidiameter" );

    previous_z = 0.0;	/* i.e.: at the entry aperture	*/

    for(
      baffle_index = num_lead_baffles - 1;
      baffle_index >= 0;
      --baffle_index )
	{
	(void)fprintf(
	  stdout,
	  "%3d %15.6f %15.6f %15.6f\n",
	  num_lead_baffles - baffle_index,
	  lead_z[ baffle_index ] - previous_z,
	  lead_z[ baffle_index ],
	  lead_y[ baffle_index ] );
	previous_z = lead_z[ baffle_index ];
	}

    for(
      baffle_index = 0;
      baffle_index < num_trail_baffles;
      ++baffle_index )
	{
	(void)fprintf(
	  stdout,
	  "%3d %15.6f %15.6f %15.6f\n",
	  num_lead_baffles + 1 + baffle_index,
	  trail_z[ baffle_index ] - previous_z,
	  trail_z[ baffle_index ],
	  trail_y[ baffle_index ] );
	previous_z = trail_z[ baffle_index ];
	}

    (void)fprintf(
      stdout,
      "    %15.6f\n",
      tube_length - previous_z );

    (void)fprintf( stdout, "\n" );

    (void)fprintf(
      stdout,
      "    Extra overlap: Entry %15.6f Exit %15.6f\n",
      0.0 - lead_shadow_z,
      trail_shadow_z - tube_length );

    exit( 0 );
    }

static double reflection_point_z(
  const double yf,
  const double y0,
  const double z0,
  const double y1,
  const double z1 )
    {
    /* Find the point at y = yf, where the angle of	*/
    /* incidence equals the angle of reflection.	*/

    double delta_y0 = yf - y0;
    double delta_y1 = yf - y1;
    double full_delta_z  = z1 - z0;
    double delta_z;

    /* Use the fact that the two right triangles are	*/
    /* similar and that the corresponding sides are	*/
    /* delta_y0 and delta_y1; combine that with the	*/
    /* fact that the sum of the two opposite legs is	*/
    /* the full_delta_z value.				*/
    /*							*/
    /*	delta_z	      full_delta_z - delta_z		*/
    /* ---------- =  ------------------------		*/
    /*  delta_y0              delta_y1			*/
    /* 							*/
    /* (Solve for delta_z.)				*/

    delta_z = full_delta_z *
      ( delta_y0 / ( delta_y0 + delta_y1 ) );

    return z0 + delta_z;
    }

static double baffle_semidiameter_y(
  const double zf,
  const double y0,
  const double z0,
  const double y1,
  const double z1 )
    {
    /* Find the y value at the baffle, given that it	*/
    /* must lie on the line connecting the two points.	*/

    double full_delta_z = z1 - z0;
    double delta_z      = zf - z0;
    double full_delta_y = y1 - y0;
    double delta_y;

    /* Just consider that at the given point, the	*/
    /* proportion of change in y must be equal to the	*/
    /* the proportion of change in z.			*/

    delta_y = full_delta_y * ( delta_z / full_delta_z );

    return y0 + delta_y;
    }

static double intersection_z(
  const double y0a,
  const double z0a,
  const double y1a,
  const double z1a,
  const double y0b,
  const double z0b,
  const double y1b,
  const double z1b )
    {
    /* Find the intersection of two lines.			*/
    /* It is assumed here that the slopes are valid.		*/

    double slope_a = ( y1a - y0a ) / ( z1a - z0a );
    double slope_b = ( y1b - y0b ) / ( z1b - z0b );
    double delta_slope = slope_a - slope_b;
    double delta_y     = y0a - y0b;

    return ( slope_a * z0a - slope_b * z0b - delta_y ) / delta_slope;
    }

static double lead_point_z(
  const double yf,
  const double y0,
  const double z0,
  const double y1,
  const double z1 )
    {
    /* Find the y-value corresponding z-value, given that it	*/
    /* must lie on the line extension of the line connecting	*/
    /* the two points.						*/

    double slope = ( y0 - y1 ) / ( z0 - z1 );
    double delta_z = ( yf - y0 ) / slope;

    return z0 + delta_z;
    }