// HDL_ANT.CPP     copyright 1994 Paul Wade N1BWT 
// version 1.0e   18 August 1994 
// e version for the metric-impaired
//
//      Any licensed radio amateur may use this program without 
//      charge; all other persons must send $73 to the 
//      ARRL Foundation, 225 Main St., Newington, CT 06111

//The source file, HDL_ANT.CPP, was compiled using Borland C++ 
//version 3.1.  Any part of it may be modified or used as part
//of another program by any amateur radio operator, provided
//that credit is given and provided there is no charge for
//its use.

//Send problems, suggestions, and comments to Paul Wade, N1BWT.
//Expect support commensurate with the cost of the program!

//antenna.h - header file for antenna.cpp
//
// N1BWT 1994
// compiled with Borland C++ version 3.1

#include <math.h>
#include <iostream.h>
#include <fstream.h>
#include <iomanip.h>
#include <ctype.h>
#include <strstrea.h>
#include <string.h>
#include <bios.h>
//#include <constrea.h>
#include <conio.h>


#include <ps.h>         // PostScript commands


#define LINESIZE 1024
#define CTRL_C 0x2e03
#define UP_ARROW 0x4800
#define DOWN_ARROW 0x5000
#define HOME_KEY 0x4700
#define END_KEY 0x4f00
#define ENTER_KEY 0x1c0d



// global constants
    const double c = 299792456.0;     //  meters per second 
    const double pi = 3.1415927 ;
    const double crad = 57.29578 ;

// other globals
    
    int lib_init = 1;		/* if lib_init, then we need to initialize
				   the library routines */
    int ph_init = 1;

    char *units = " mm.  ";

    float metric = 1.0;

    char *unit_type = "Metric ";
//






// start of ANTENNA.H - class declarations

class Antenna           // Abstract base class
{
public:
    Antenna();
    void helpscreen();
    void get_freq();
    void get_wg();
    void horn_dimensions();
    void horn_calcgain();
    void horn_phasecenters();
    void horn_angles();
    void banner(ostream&);
    void info_listing(ostream&);
    void post_print(ostream&);
    void PANFI_corr();
    void range_calc();
    int bad_input() {return bad;};
      
protected:
    double freq;      
    double H_wg;                // waveguide H dim. in mm. */ 
    double E_wg;                // waveguide E dim. in mm. */ 
    double wavelength;          // in mm */
    double axial_length;        // axial length of feed horn flare in mm */
    double H_aperture;          // horn E-Plane aperture in mm */ 
    double E_aperture;          // horn H-Plane aperture in mm */ 
    double horn_gain;       
    double E_phase_center;
    double H_phase_center;
    double E_flare_angle;       // angle from axis to E-side
    double H_flare_angle;       // angle from axis to H-side
    double E_side;
    double E_ang;
    double H_side;
    double H_ang;
    double aperture;
    double efficiency;
    double rayleigh_dist;
    double directivity;
    int bad;                    // problem with input
};


class Horn : public Antenna
{
public:
    Horn();
    void lookup();
    void simple_gain();
    void coords();
    void PS_template();
    void write_data(ostream&, int);
        
private:

    double simplegain;       
    Point inner[6];
    Point outer[6];

};
//







class Dish : public Antenna
{
public:
    Dish();
    void d_init();
    void illumination();
    void calc_gain();
    void phasecenters();
    void PS_template();
    
private:
                // all dimensions in mm.
    double freq;
    double diameter;                
    double depth;                   // at center of dish
    double focal_length;        
    double f_over_D;
    double feedangle;               // diameter subtended angle from focal point
    double feedang_deg;
    double R_edge;                  // distance from focal point to edge
    double space_attn;
    double desired_taper;
    double BW3;                     // 3 dB beamwidth
    double simple_horn_diameter;
    double dish_gain;                 
};

class Lens : public Antenna
{

public:
    Lens();
    void dimension();    
    void calculate();
    void phasecheck();
    void plates();
    void write_data(ostream &, int);
    void crude_graphics(ostream &);
    
private:

    double Lens_diam ;       /* Lens diameter in mm */
    double spacing   ;       /* metal plate spacing in mm */
    double width_increase ;  /* plates get wider toward edge */
    double shift[25];        /* H-plane phase compensation */
    int N ;                  /* number of plates in lens */
    double BW_E_horn ;           /* feed horn E-plane Beamwidth in degrees */
    double lens_focal_length ;   /* in mm */
    double Gain_lens ;           /* dB added by lens */
    double index ;               /* effective index of refraction */
    double lens_to_horn;        // horn mouth to lens curve center in mm.
    double step;                // step distance in zoned lens plate */
    int compensate ;            // flag for phase center compensation needed
};

// end of ANTENNA.H - class declarations

//






// start of ANTENNA.CPP - class member functions

  
Antenna::Antenna()              // Constructor
{
    freq = 0.0;      
    H_wg = 0.0; 
    E_wg = 0.0;
    wavelength = 0.0;
    freq = 0.0;        
    H_wg = 0.0;        
    E_wg = 0.0;        
    wavelength = 0.0;  
    horn_gain = 0.0;        
    axial_length = 0.0;
    H_aperture = 0.0;
    E_aperture = 0.0;   
    H_flare_angle = 0.0;
    E_flare_angle = 0.0;
    E_side = 0.0;
    E_ang = 0.0;
    H_side = 0.0;
    H_ang = 0.0;
    aperture = 0.0;
    efficiency = 0.0;
    rayleigh_dist = 0.0;
    directivity = 0.0;
    bad = 0;
}

void Antenna::banner(ostream& ostrm)
{  /* a banner */
        
    ostrm << "\n\n\n\n";
    ostrm << "             **************** \n";
    ostrm << "             *              * \n";
    ostrm << "             *    HDL_ANT   * \n";
    ostrm << "             *              * \n";
    ostrm << "             **************** \n";
    ostrm << "\n\n     copyright 1994 Paul Wade N1BWT \n\n";
    ostrm << "      Version 1.0E with option for the metric-impaired.\n\n";
    
    ostrm << "      Any licensed radio amateur may use this program "
          << "without \n";
    ostrm << "      charge; all other persons must send $73 to the \n";
    ostrm << "      ARRL Foundation, 225 Main St., Newington, CT 06111";
    ostrm << "\n\n";
}

void Antenna::info_listing(ostream& ostrm)
{

ostrm << "HDL_ANT does the calculations needed to successfully implement\n";
ostrm << "several microwave antennas: horns, parabolic dishes, and \n";
ostrm << "metal-plate lenses.  It also does the calculations involved in\n";
ostrm << "setting up an antenna range for antenna gain measurement.  More\n";
ostrm << "detailed explanations of the antennas and measurements may be\n";
ostrm << "found in the series of QEX articles entitled 'Practical Microwave\n";
ostrm << "Antennas,' in the September, October, and November 1994 issues.\n\n";

ostrm << "All dimensions are in millimeters, since one mm. is about the \n";
ostrm << "accuracy required for antennas to work well at 10 GHz.\n";
ostrm << "However, there is an option for dimensions in inches, for the\n";
ostrm << "metric-impaired.\n\n";

ostrm << "MENU CHOICES:\n\n";
ostrm << "  H\n\n";

ostrm << "For horn antennas, HDL_ANT will design an 'optimum' horn with \n";
ostrm << "a specified gain (between 10 and 25 dB) for any frequency - use \n";
ostrm << "menu entry 'H' for this option.  It will then create a PostScript\n";
ostrm << "template file which may be used to print a paper template for\n";
ostrm << "horn construction.  The paper template is attached to a sheet of\n";
ostrm << "copper or brass which is cut and folded on the template lines,\n";
ostrm << "then soldered together and soldered to a waveguide.\n\n";
ostrm << "See menu entry 'P' for more details on printing a "
      << "PostScript file.\n\n";

ostrm << "  E\n\n";

ostrm << "Menu entry 'E' is used to calculate the gain and other parameters\n";
ostrm << "for an existing horn - perhaps one you found at a hamfest.  It \n";
ostrm << "can also create a template, to duplicate a horn known to have\n";
ostrm << "good dimensions.\n\n";

ostrm << "  D\n\n";

ostrm << "Parabolic dish calculations are under menu entry 'D'.  The most\n";
ostrm << "important dish dimension is the focal length, and the easiest way\n";
ostrm << "to find it is to measure the depth of the dish in the center,\n";
ostrm << "measured from a straightedge across the rim of the dish.  HDL_ANT\n";
ostrm << "then calculates the focal length and f/D ratio, to allow\n";
ostrm << "selection of the best feed antenna - see the QEX series for\n";
ostrm << "details on feedhorn selection.  For smaller dishes, HDL_ANT can\n";
ostrm << "also generate a PostScript template file;  the printed template\n";
ostrm << "is used to check the accuracy of the parabolic curve of the dish.\n";
ostrm << "If it is not a good fit, sometimes a template for a slightly\n";
ostrm << "larger or smaller f/D ratio gives a better fit, since dishes\n";
ostrm << "sometimes have a flat area or large hole in the center which\n";
ostrm << "yields a misleading depth measurement.\n\n";

ostrm << "  L\n\n";

ostrm << "Menu entry 'L' is used to design metal-plate lens antennas.\n";
ostrm << "These are normally fed by a horn; if an appropriate horn is not\n";
ostrm << "already available, then the 'H' entry may be used to design the\n";
ostrm << "horn.  It may take a few trials to find a good combination of\n";
ostrm << "horn and lens.  Note that the 'optimum' horns designed under the\n";
ostrm << "'H' entry of the HDL_ANT program do not usually have the matched\n";
ostrm << "E-plane and H-plane phase centers needed for good lens\n";
ostrm << "performance; some adjustment of horn dimensions can improve the\n";
ostrm << "phase center matching.  The lens design does not create templates\n";
ostrm << "since the lens curves are all simple circles which can be drawn\n";
ostrm << "more accurately with a compass.\n\n";

ostrm << "  R\n\n";

ostrm << "Antenna range design is selected by menu entry 'R'.  HDL_ANT\n";
ostrm << "calculates the minimum length and measurement height required for\n";
ostrm << "a specified maximum antenna aperture, and then calculates the\n";
ostrm << "height of the source antenna needed to have the RF field peak at\n";
ostrm << "the measurement height.\n\n";

ostrm << "  C\n\n";

ostrm << "If a noise figure meter is used to measure antenna gain as\n";
ostrm << "described in the QEX articles, then it is necessary to convert\n";
ostrm << "noise figure readings into antenna gain.  Menu entry 'C' makes\n";
ostrm << "these calculations.\n\n";

ostrm << "  P\n\n";

ostrm << "Menu entry 'P' is a short explanation of how to print PostScript\n";
ostrm << "files.\n\n";

ostrm << "  I\n\n";

ostrm << "Menu entry 'I' is this summary of information about HDL_ANT.\n\n";

ostrm << "  U\n\n";

ostrm << "Menu entry 'U' allows you to select dimensions in either\n" 
      << "millimeters or inches.  HDL_ANT defaults to millimeters.\n\n";

ostrm << "  Q\n\n";

ostrm << "'Q' for Quit ends the program.\n\n\n";


ostrm << "The source file, HDL_ANT.CPP, was compiled using Borland C++ \n";
ostrm << "version 3.1.  Any part of it may be modified or used as part\n";
ostrm << "of another program by any amateur radio operator, provided\n";
ostrm << "that credit is given and provided there is no charge for\n";
ostrm << "its use.\n\n";

ostrm << "Send problems, suggestions, and comments to Paul Wade, N1BWT.\n";
ostrm << "Expect support commensurate with the cost of the program!\n";
}

void Antenna::post_print(ostream& ostrm)
{
ostrm << "\n            PRINTING POSTSCRIPT[tm] FILES\n\n";
        
ostrm << "      The easiest way to print PostScript files is with a \n";
ostrm << "      PostScript compatible laser printer.  These have become\n";
ostrm << "      more affordable and are becoming more common; for\n";
ostrm << "      instance, the public library in my small town has one\n";
ostrm << "      attached to a public-access computer.   However, they are\n";
ostrm << "      still roughly twice as expensive as the dot-matrix\n";
ostrm << "      printers that most of us use with our personal computers.\n";
ostrm << "      \n";
ostrm << "      An alternative to a laser printer is software that\n";
ostrm << "      interprets PostScript language commands for display on a\n";
ostrm << "      computer VGA display or a dot-matrix printer.  I know of\n";
ostrm << "      several versions of this type of software:  GoScript, a\n";
ostrm << "      commercial product, Ultrascript, a commercial product,\n";
ostrm << "      Freedom of Press , a commercial product, and Ghostscript,\n";
ostrm << "      [Copyright (C) 1990, 1992 Aladdin Enterprises.  All\n";
ostrm << "      rights reserved.  Distributed by Free Software\n";
ostrm << "      Foundation, Inc.], which is freeware.  \n";
ostrm << "      \n";
ostrm << "      I have only used Ghostscript, version 2.5.  While 286,\n";
ostrm << "      386, and windows versions are available, the 286 version\n";
ostrm << "      seems to work most reliably even on a 386 or 486 PC.  It\n";
ostrm << "      uses Unix style command strings which are difficult to\n";
ostrm << "      remember, so I've included two BAT files to help:\n";
ostrm << "      GS_VIEW.BAT for viewing on a screen, and GS_PRINT.BAT for\n";
ostrm << "      printing on an Epson dot-matrix printers.  For other\n";
ostrm << "      brands of printer, the command will have to be changed\n";
ostrm << "      appropriately, which will require reading of the\n";
ostrm << "      documentation.  Type GS_VIEW <filename.ps> or GS_PRINT\n";
ostrm << "      <filename.ps> to use them.  Be sure to type QUIT when you\n";
ostrm << "      are through or your PC may be left in an unhappy state\n";
ostrm << "      requiring rebooting.  \n";
ostrm << "      \n";
ostrm << "      I've also included a sample PostScript file, SQUARE.PS,\n";
ostrm << "      which draws a four inch square to make sure that\n";
ostrm << "      templates will be drawn to scale, and a sample horn\n";
ostrm << "      template, HORN18.PS, to get you started.  If the\n";
ostrm << "      dimensions of the printed square are slightly off, you\n";
ostrm << "      can correct them.  Each template has a line near the\n";
ostrm << "      beginning of the file:\n";
ostrm << "      \n";
ostrm << "              1.0 1.0 scale\n";
ostrm << "              \n";
ostrm << "      the first number is the scale factor in the x\n";
ostrm << "      (horizontal) direction, and the second is the scale\n";
ostrm << "      factor in the y (vertical) direction.  Edit the SQUARE.PS\n";
ostrm << "      file with an editor to change these numbers slightly;\n";
ostrm << "      when you find a combination that prints a square exactly\n";
ostrm << "      four inches on a side, then you have compensated for your\n";
ostrm << "      printer.  Edit these same numbers into any template to be\n";
ostrm << "      printed on the same printer and the dimensions will come\n";
ostrm << "      out right.\n";
ostrm << "      \n";
ostrm << "      The Ghostscript files are available on many bulletin\n";
ostrm << "      boards and Internet locations.  They are in ZIP format,\n";
ostrm << "      so they must me downloaded, unZIPped, and installed\n";
ostrm << "      according to the README documentation.\n";
ostrm << "      \n";
ostrm << "      I have not used any of the commercial products, but I\n";
ostrm << "      would expect a commercial product to be much easier to\n";
ostrm << "      install and use than freeware or shareware.   \n";
ostrm << "      \n";
ostrm << "      The batch files are as follows:\n\n";
ostrm << "GS_VIEW.BAT\n\n";
ostrm << "   gs %1\n\n";
ostrm << "GS_PRINT.BAT\n\n";
ostrm << "   gs -sDEVICE=epson -r60x60 %1\n";
}

void Antenna::PANFI_corr()
{
// correction for PANFI readings to Antenna Gain

    float Tex_dB;  // Excess noise calibration

    cout << "\nPANFI calibration setting - Excess noise in dB: " ;
    cout.flush();
    cin >> Tex_dB;

    float std_db;
    cout << "\nGain of standard gain antenna: " ;
    cout.flush();
    cin >> std_db;
    cout << "\n\n";

    double nf = 0.0;   // indicated noise figure in dB

    cout << "Enter a negative Noise Figure to return to Main Menu\n\n\n";

    while (nf >= 0.0)
    {
        double std_nf;
        cout << "\nIndicated Noise Figure for Standard gain antenna: " ;
	cout.flush();
	cin >> std_nf ;

        double ratio = pow(10.0,(0.1 * (Tex_dB - std_nf)));
        double Y_std = 10.0 * log10( 1.0 + ratio);
        cout << setprecision(2) 
	     << "\n            Y factor = " << Y_std <<  " dB\n\n";

        cout << "\nIndicated Noise Figure for antenna under test: " ;
	cout.flush();
	cin >> nf ;

        ratio = pow(10.0,(0.1 * (Tex_dB - nf)));
        double Y_db = 10.0 * log10( 1.0 + ratio);

        cout << setprecision(2) 
	     << "\n           Y factor = " << Y_db <<  " dB\n\n";


	double G_db = Y_db - Y_std + std_db;
	cout << setprecision(2)
	     << "      Corrected gain = " << G_db <<  " dBi\n\n";

	float dia;
	cout << "Aperture diameter in " << units << ": " ;
	cout.flush();
	cin >> dia ;
        dia *= metric;

        double g_ratio = pow(10.0,(0.1 * G_db));

        directivity = pi * pi * dia * dia /( wavelength * wavelength);

//cout << "Directivity = " << directivity;

        efficiency = g_ratio / directivity;
                
        cout << setprecision(2) 
             << "\n\n      Aperture efficiency = " << (efficiency * 100) 
             <<  " %\n\n\n";        
    }
}


void Antenna::range_calc()
{
// calculations for Antenna Range design
// this routine returns feet or meters so has hard constants

    cout << "\nLargest aperture diameter to be measured, in "
         << units << ": ";
    cout.flush();
    cin >> aperture;
    
    aperture *= metric;

    rayleigh_dist = 2.0 * aperture * aperture / wavelength;

    float big_u;
    char *big_units = "xxxxxx";

    if (metric == 1.0)
    {
        big_u = 1000.0;
        big_units = "meters";
    }
    else
    {
        big_u = 304.8;
        big_units = "feet  ";
    }
        
     cout << "\n\nRayleigh distance (minimum range length) = "
          << setprecision(2) << (rayleigh_dist / big_u) << " "
          << big_units << "\n\n";

     cout << "Suggested measurement height = "
          << setprecision(2) << (aperture * 4.0 / big_u) << " "
          << big_units << "\n\n";

    double range_length;
    cout << "Enter desired range length in " << big_units << ": ";
    cout.flush();
    cin  >> range_length;
    
    range_length *= big_u;      // mm.
    
    if (range_length < rayleigh_dist)
        cout << "\n\nWARNING - range too short for "
             << "accurate gain measurement\n\n";

        double meas_height;
        cout << "\n\nEnter desired measurement height in "
             << big_units << ": ";
        cout.flush();
        cin  >> meas_height;
             
        meas_height *= big_u;      // mm.
    
    if (meas_height < (3.2 * aperture))
        cout << "\n\nWARNING - measurement height insufficient for "
             << "accurate gain measurement\n\n";

    double source_height = wavelength * range_length / (4.0 * meas_height);
    
    cout << "\n\nHeight of source antenna should be "
         << setprecision(2) << (source_height / metric) << " "
         << units << " \n\n";
}


//




Dish::Dish() : Antenna()
{               // constructor 
    freq = 0.0;                 
    diameter = 0.0;             
    depth = 0.0;                
    focal_length = 0.0;         
    f_over_D = 0.0;             
    feedangle = 0.0;            
    feedang_deg = 0.0;          
    R_edge = 0.0;                    
    space_attn = 0.0;           
    desired_taper = 0.0;        
    BW3 = 0.0;                  
    simple_horn_diameter = 0.0; 
    wavelength = 0.0;           
    dish_gain = 0.0;                 
    
}

void Dish::d_init()
{
    cout << "\nDiameter of dish in " << units << ": " ;
    cout.flush();
    cin >> diameter ;
    diameter *= metric;   //mm

    char mode = 'n';
    cout << "\nDo you know the f/D ratio [Yes or No]?  ";
    cout.flush();
    cin >> mode ;
    mode = toupper(mode);

    if (mode == 'Y')
    {
        cout << " \n Enter f/D ratio: ";
	cout.flush();
	cin >> f_over_D;

	focal_length = f_over_D * diameter ;
    }
    else
    {
	cout << " \n Enter Depth of dish in " << units << ": " ;
	cout.flush();
	cin >> depth;
        depth *= metric;      //mm

	focal_length = (diameter * diameter) / ( 16.0 * depth);
	f_over_D = focal_length / diameter ;

	cout << "\n\n     f/D = " << setprecision(2) << f_over_D << "\n" ;
    }

    cout << "\nFocal length = "
         << setprecision(2) <<  (focal_length / metric) 
         << " " << units << endl;

}

void Dish::illumination()
{
	feedangle = 2 * atan2(1.0,((2.0 * f_over_D) - (1.0/(8.0 * f_over_D))));
	feedang_deg = crad * feedangle ;   // degrees conversion

		// R_edge = distance from focal point to edge

	R_edge = 2.0 * focal_length / ( 1 + cos(feedangle/2.0));

		// Space_attn is additional path loss to edge of dish
	double Space_attn = (20.0 * log10(R_edge / focal_length));

		// Desired_taper to achieve 10 dB actual edge taper
	double Desired_taper = 10.0 - Space_attn;  // dB

	BW3 = feedangle * sqrt( 3.0 / Desired_taper) ;

	double Simple_horn_diameter = 66.0 / (BW3 * crad) ; // wavelengths

	cout << "\n\n  Illumination angle for feed = "
	     << setprecision(2) <<  feedang_deg << " degrees";

	cout << "\n  Space attenuation = "
	     << setprecision(2) << Space_attn << " dB ";
	cout << "\n  Desired taper = "
	     << setprecision(2) << Desired_taper << " db "
	     << "for 10 dB edge illumination";

	cout << "\n\nA simple feed horn would have a diameter of "
	     << setprecision(2) << Simple_horn_diameter << " wavelengths";

	cout << "\n   for a 3 dB beamwidth of "
	     << setprecision(2) << (BW3 * crad) << " degrees\n";


//        Write_PS(diameter, focal_length);

}

void Dish::calc_gain()
{
    efficiency = 0.5;     // assume 50% efficiency

    double radius = 0.5 * diameter;

    aperture = pi * radius * radius;

    dish_gain = efficiency * aperture * 4.0 * pi / (wavelength * wavelength);

    dish_gain = 10.0 * log10(dish_gain);

    cout << "\nGain at 50% efficiency = " << setprecision(2) << dish_gain
	 << " dBi \n";
}


void Dish::PS_template()
{                      // PostScript template on one page

    char descriptor[LINESIZE] = " ";
    char outname[LINESIZE];

    cin.ignore(LINESIZE,'\n');
    cout << " Dish description line: ";
    cout.flush();
    cin.get(descriptor,LINESIZE,'\n');

    cout << "File name for dish template (should be .PS ): ";
    cout.flush();
    cin >> outname ;

    ofstream outfil(outname);

    if (!outfil)
    {
	cerr << "Cannot open output file for output";
    }

// header

    outfil << "%! PostScript\n";

    outfil << "% template for " << (diameter/metric) 
           << " " << units << " parabolic dish\n";
    outfil << "% with f/d = " << f_over_D << " and focal length = "
           << (focal_length/metric) << " " << units << "\n";
    
    outfil << "% generated by HDL_ANT by N1BWT 1994 \n";
    outfil << "%\n";
    outfil << "%" << descriptor;
    outfil << "%\n";
    outfil << "% Print on PostScript printer\n";
    outfil << "% cut along the curve,\n";
    outfil << "% and use it as a template for the dish.\n\n\n";
    
    outfil << "%\n";
    outfil << "/in { 72 mul } def\n";
    outfil << "/mm { 2.834646 mul } def\n";

    outfil << "%\n";
    outfil << "% 100% scaling factors - change for printer correction\n";
    outfil << "1.0 1.0 scale \n";
    outfil << "%\n";

outfil << "% DISH **********************************************\n\n";

//    Point offset =  Point(50.0,(100.0 - (diameter/4.0))) ;
    Point offset =  Point(50.0, 260.0);

    outfil << offset.x() << " mm " << offset.y() << " mm translate" << endl;

    outfil << "% DISH\n";
outfil << "gsave\n";

// and dish outline 

// NOTE: I'm doing this directly in picas, so can't use PS_lineto
        
    int prad = int(2.834636 * diameter / 2.0);         // radius in picas 
    double pa = 4.0 * focal_length * 2.834636;       // 4 * f in picas 
        
    outfil << "newpath\n";

    outfil << PS_moveto() << endl;
        
    double i_d;    
    double curvature;
    for ( int i = 0; i <= prad; i++)
    {
        i_d = (double)i;
        curvature = i_d * i_d / pa ;
        outfil << curvature << " " << (-i) << " lineto\n";   // draw it vertically
    }
    outfil << "stroke \n\n";
            
// and draw a box around it 

    outfil << "newpath\n";
    outfil << curvature << " " << (-prad) << " moveto" << endl;
    outfil << curvature << " 0 lineto " << endl;
    outfil << "0 0 lineto " << endl;
    outfil << "0 " << (-prad) << " lineto" << endl;
    outfil << curvature << " " << (-prad) << " lineto" << endl;
    outfil << "stroke\n";    

// Text
    outfil << "/Helvetica-Bold findfont\n";
    outfil << "14 scalefont\n";
    outfil << "setfont\n";
// position title

    outfil << "-25 mm -25 mm moveto\n";
    outfil << "(Template for " << diameter << "mm dish with f/d = " 
           << f_over_D << " ) show\n";
    outfil << "-25 mm  -35 mm moveto\n";
    outfil << " (N1BWT 1994) show \n";


    outfil << "% END DISH **********************************************\n";
    outfil << "grestore\n";

    outfil << "showpage\n";
    outfil << "%\n";

    outfil.close();
//cerr << "outfile closed";
}
 


//






Horn::Horn() : Antenna()
{               // constructor 
    freq = 0;        
    H_wg = 0;        
    E_wg = 0;        
    wavelength = 0;  
    horn_gain = 0;        
    axial_length = 0;
    H_aperture = 0;
    E_aperture = 0;   
    E_phase_center = 0;
    H_phase_center = 0;
    simplegain = 0;        
    E_flare_angle = 0;
    E_side = 0;
    E_ang = 0;
    H_flare_angle = 0;
    H_side = 0;
    H_ang = 0;
}
            

void Antenna::get_freq()
{                   // reads initial data
    cout << "\n" << "Frequency in MHz: " ;
    cout.flush();
    cin >> freq ;

    wavelength = c / (1000.0 * freq);
}

void Antenna::get_wg()
{
    cout << "         ____________     \n";
    cout << "        |     E      |   waveguide  \n";
    cout << "        |     ^ -->H |             Hint: for WR-90\n";
    cout << "        |     |      |               H = 22.86 mm. = 0.9 in.,\n";
    cout << "         ------------                E = 10.16 mm. = 0.4 in. \n\n";

    bad = 1 ;
    while (bad >= 1)
    {    
        cout << "\n" << "Enter H-plane inside dimension of waveguide in "
             << units << ": ";
	cout.flush();
	cin >> H_wg;
        H_wg *= metric;         //mm

        cout << "\n" << "Enter E-plane inside dimension of waveguide in "
             << units << ": ";
	cout.flush();
	cin >> E_wg;
        E_wg *= metric;         //mm
    
	if (wavelength < (0.98 * H_wg) || wavelength > ( 1.72 * H_wg ))
                        /* slight fudge to get edges */
        {
            cout << " ERROR - frequency out of waveguide range \n";
            bad++ ;
            if (bad > 3) break;
        }
        else bad =0;
    }
}

void Antenna::horn_dimensions()
{                               // input dimensions of a horn
        
    cout << "\n" << "ENTER DESIRED PHYSICAL DIMENSIONS OF HORN: ";
    

    cout << "\n" << "Enter axial length of horn in "
         << units << ": ";
    cout.flush();
    cin >> axial_length;
    axial_length *= metric;
    
    cout << "\n" << "Enter H-plane aperture of horn in "
         << units << ": ";
    cout.flush();
    cin >> H_aperture;
    H_aperture *= metric;
    
    cout << "\n" << "Enter E-plane aperture of horn in "
         << units << ": ";
    cout.flush();
    cin >> E_aperture;
    E_aperture *= metric;
}        

void Horn::lookup()
{

//   Lookup table based on 
//    D.E. Cozzens, "Tables Ease Horn Design," _Microwaves_, March 1966, p37.
//   from which I eyeballed nominal values;
//   anyway, the program calculates a new gain later 
         
    float G[ 26];
    
    for (int i = 0; i < 10; i++)
        G[i] = 0;
        
    G[10] = 0.080;
    G[11] = 0.24;
    G[12] = 0.4;
    G[13] = 0.6;
    G[14] = 0.88;
    G[15] = 1.25;
    G[16] = 1.7;
    G[17] = 2.25;
    G[18] = 3.0;
    G[19] = 4.0;
    G[20] = 5.2;
    G[21] = 6.7;
    G[22] = 8.7;
    G[23] = 11.2;
    G[24] = 14.3;
    G[25] = 18.3;

    bad = 1;
    
    double gain;    

    while (bad)
    {
        cout << "\n" << "Enter desired gain of horn in dB: ";
	cout.flush();
	cin >> gain;

        if ((gain < 10.0) || (gain > 25.0))
            cout << "Gain must be between 10 and 25 dB";
        else bad = 0;
    }        

    axial_length = wavelength * G[int(gain)];

// {need gain as a ratio to find horn apertures }
// {**** NOTE y=z**x -> y:=exp(x*ln(z)) ****}

    double gain_ratio = exp(0.2302585 * (int(gain)));            

    H_aperture = wavelength * 0.4675 * sqrt(gain_ratio);

    E_aperture = wavelength * 0.3463 * sqrt(gain_ratio);

    cout << "\n" << " CALCULATED HORN DIMENSIONS: \n";   
    
    cout << "\n" << "   Axial length of horn is "
         << setprecision(2) << (axial_length/metric) << " "
	 << units;
    cout << "\n" << "   H-plane aperture of horn is "
	 << setprecision(2) << (H_aperture/metric) << " "
	 << units;
    cout << "\n" << "   E-plane aperture of horn is "
	 << setprecision(2) << (E_aperture/metric) << " "
	 << units;
}


void Horn::simple_gain()
{                 // temp - simple formula from Kraus

    simplegain = 4.5 * (H_aperture/wavelength) * (E_aperture/wavelength);

    simplegain = 10.0 * log10(simplegain) + 2.14 ;   //{ dBi }
            
//    cout << " Simple gain of horn is "
//         << setprecision(2) << simplegain << " dBi \n";
}            

double Le(double s)
{
//    MathCad fit by Matt Reilly of Fig. 23 from Balanis

//    input is Maximum input phase error in wavelengths
//    returns loss due to phase error in dB.

//	Le(s) = sum(from i=1 to 10) Xle[i] * s^i

    if ( s < 0 || s > 1)
    {
        cerr << "ERROR: excessive phase error in Le\n\n";
        return 1;
    }

    double Xle[11];

// Xle[i] = 
    Xle[1] =   10.5075396618;
    Xle[2] =  -236.610713851;
    Xle[3] =   2603.14290758;
    Xle[4] =  -15318.8133545;
    Xle[5] =   55318.7209697;
    Xle[6] =  -125471.932662;
    Xle[7] =   178426.752358;
    Xle[8] =  -153811.177005;
    Xle[9] =   73233.5757222;
    Xle[10] =  -14743.1657614;

    double Le = 0;

    for (int i = 1; i < 11; i++)
    {
	Le += Xle[i] * pow(s,i);
    }
    return Le;
}

double Lh(double t)
{
//    MathCad fit by Matt Reilly of Fig. 23 from Balanis

//    input is Maximum input phase error in wavelengths
//    returns loss due to phase error in dB.

//	Lh(t) = sum(from i=1 to 10) Xlh[i] * t^i

    if ( t < 0 || t > 1)
    {
        cerr << "ERROR: excessive phase error in Lh\n\n";
        return 1;
    }

    double Xlh[11];

// Xlh[i] =
    Xlh[1] =   -13.9134920465;
    Xlh[2] =      398.4347217;
    Xlh[3] =   -4241.37951372;
    Xlh[4] =    24010.7335428;
    Xlh[5] =    -79636.863311;
    Xlh[6] =    162957.754379;
    Xlh[7] =   -208325.065791;
    Xlh[8] =    162037.036743;
    Xlh[9] =   -70133.3772858;
    Xlh[10] =     12951.940007;

    double Lh = 0;

    for (int i = 1; i < 11; i++)
    {
        Lh += Xlh[i] * pow(t,i);
    }
    return Lh;
}

double phase_error(double aperture, double flare_angle, double wavelength)
{
    double rho = 0.5 * aperture / sin(flare_angle);
//cout << setprecision(4) << "\n aperture = " << aperture
//     << " flare angle = " << flare_angle
//     << " rho =  " << rho ;   
    double error = aperture * aperture / ( 8.0 * wavelength * rho );
    
    return error;
}

void Antenna::horn_calcgain()
{
//improved algorithm from
// Balanis, C.A., "Horn Antennas", Chapter 8 in 
// Lo, Y.T. and Lee, S.W.
// Antenna Handbook: Theory, Applications and Design, Van Nostrand, 1988
// pp 8-39 to 8-42.


    horn_gain = 10.0 * ( 1.008 + 
                 log10( E_aperture * H_aperture / (wavelength * wavelength)));


// correct for phase error


    double s = phase_error(E_aperture, H_flare_angle, wavelength);
//cout << setprecision(4) << "\n s = " << s << "  Le(s) = " << (Le(s));

    double t = phase_error(H_aperture, E_flare_angle, wavelength);
//cout << setprecision(4) << "\n t = " << t << "  Lh(t) = " << (Lh(t)) << endl;

    horn_gain -= (Le(s) + Lh(t));               // dBi

//    cout << " Improved gain of horn is "
//         << setprecision(2) << horn_gain << " dBi \n";

}


void Horn::write_data(ostream& ostrm, int interactive)
{
    // anti warning hack.
    int i, j;
    for(i = j = 0; i < 10; i++) j = j + interactive;

    ostrm << " \nSimple gain of horn is "
	  << setprecision(2) << simplegain << " dBi \n\n";

    ostrm << " Improved gain of horn is "
	  << setprecision(2) << horn_gain << " dBi \n\n";

    ostrm << "H-Plane phase center is "
         << setprecision(2) << H_phase_center
         << " wavelengths inside horn mouth\n\n"; 
    ostrm << "E-Plane phase center is "
	  << setprecision(2) << E_phase_center
	  << " wavelengths inside horn mouth\n";
}

void Antenna::horn_angles()
{                
// note definition of E_flare_angle is the angle that the E side is 
// flared at, which is the flare angle in the H plane

    E_flare_angle = atan2((0.5 * (H_aperture - H_wg)) , axial_length);
    E_side = axial_length / cos(E_flare_angle);
    E_ang = atan2((0.5 * (E_aperture - E_wg)) , E_side);

    H_flare_angle = atan2((0.5 * (E_aperture - E_wg)) , axial_length);
    H_side = axial_length / cos(H_flare_angle);
    H_ang = atan2((0.5 * (H_aperture - H_wg)) , H_side);

// test

//    cout << setprecision(4) 
//         << "Flare angles " << E_flare_angle << "   " << H_flare_angle << "\n";
//    cout << setprecision(4) 
//         <<  "Side lengths" << E_side << "   " << H_side << "\n";
//    cout << setprecision(4) 
//         << "Angles " << E_ang << "   " << H_ang << "\n";
}


void Horn::coords()
{
// TEST
//    cout << " ENTERING fixed Horn::coords()" << endl;

    double fold_ang = E_ang + H_ang;

    inner[1] = inner[0] + Point(E_wg,0);

    inner[2] = inner[1] + Point((H_wg * cos(fold_ang)),-(H_wg * sin(fold_ang)));

    inner[3] = inner[2] + Point((E_wg * cos(2 * fold_ang)),
                                -(E_wg * sin(2 * fold_ang)));
                           
    inner[4] = inner[3] + Point((H_wg * cos(3 * fold_ang)),
                                -(H_wg * sin(3 * fold_ang)));

    outer[0] += Point(-(0.5 * (E_aperture - E_wg)),E_side);

    outer[1] = outer[0] + Point(E_aperture,0);

    outer[2] = outer[1] + Point((H_aperture * cos(fold_ang)),
				-(H_aperture * sin(fold_ang)));

    outer[3] = outer[2] + Point((E_aperture * cos(2 * fold_ang)),
				-(E_aperture * sin(2 * fold_ang)));

    outer[4] = outer[3] + Point((H_aperture * cos(3 * fold_ang)),
				-(H_aperture * sin(3 * fold_ang)));

	// and a solder flap


    inner[5] = inner[4] + Point((7 * cos(3*fold_ang - 0.7854 + H_ang)),
			       (-(7 * sin(3*fold_ang - 0.7854 + H_ang))));

    outer[5] = outer[4] + Point((7 * cos(3*fold_ang + 0.7854 + H_ang)),
                                (-(7 * sin(3*fold_ang + 0.7854 + H_ang))));



//%t   cout << "H_angle = " << H_ang << endl;
    
//%t   cout << "fold_angle = " << fold_ang << endl;
    
//%t   Point in5_pt = Point((7 * cos(3*fold_ang - 0.7854 + H_ang)),
//			       (-(7 * cos(3*fold_ang - 0.7854 + H_ang))));

//%t   cout << "inner_5 increment = " << in5_pt << endl;

//%t   Point out5_pt = Point((7 * cos(3*fold_ang + 0.7854 + H_ang)),
//				(-(7 * cos(3*fold_ang + 0.7854 + H_ang))));

//%t   cout << "outer_5 increment = " << out5_pt << endl;

}
//






void Horn::PS_template()
{                      // PostScript template on one page
                       // or two pages for larger horns
                       
    char descriptor[LINESIZE] = " ";
    char outname[LINESIZE];

    cin.ignore(LINESIZE,'\n');
    cout << " Horn description line: ";
    cout.flush();
    cin.get(descriptor,LINESIZE,'\n');

    cout << "File name for horn template (should be .PS ): ";
    cout.flush();
    cin >> outname ;

    ofstream outfil(outname);

    if (!outfil) 
    {
        cerr << "Cannot open output file for output";
    }
    
    int pages = 2;
    
    Point offset;

    if ((outer[0].y() - outer[4].y()) < 260)         // horn fits on page 
    {
        pages = 1;
        if (outer[4].x() < outer[0].x())
            offset = Point(-outer[4].x(), -outer[0].y());
        else offset = Point(-outer[0].x(), -outer[0].y());
    }
    else offset = Point(-outer[0].x(), -outer[0].y());

    offset =  Point(25,254) + offset;
            
    
//      PostScript header 

    outfil << "%! PostScript \n";
    outfil << "/mm { 2.834646 mul } def\n";
    outfil << "/in { 72 mul } def\n";

    outfil << "%\n";
    outfil << "% 100% scaling factors - change for printer correction\n";
    outfil << "1.0 1.0 scale \n";
    outfil << "%\n";

    outfil << setprecision(2);
    outfil << setprecision(2)
           << "% Template for " << horn_gain << " dBi pyramidal horn for" 
           << freq << " MHz" << endl;
    outfil << "% generated by HDL_ANT by N1BWT 1994 " << endl;
    outfil << "%" << endl;
    outfil << "%" << descriptor;
    outfil << "%" << endl;
    outfil << "% Print on PostScript printer" << endl;
    outfil << "% stick to flashing copper";
    outfil << "% and cut out pieces." << endl;
    outfil << "% Then fold up and solder horn";
    outfil << "% then solder it to waveguide flange." << endl;
    outfil << "% Waveguide size is " << (H_wg/metric) << " "
           << units << " x " << (E_wg/metric) << " " << units << "\n";
    outfil << "%" << endl;
    outfil << "%" << endl;
    
    outfil << "gsave" << endl;
    outfil << "% HORN **********************************************" << endl;
    outfil << "%" << endl;
    
    outfil << offset.x() << " mm " << offset.y() << " mm translate" << endl;
    
    outfil << "%" << endl;
    outfil << "%" << endl;
    outfil << "% HORN" << endl;
    outfil << "%" << endl;

// Text

    outfil << "/Helvetica-Bold findfont" << endl;
    outfil << "14 scalefont" << endl;
    outfil << "setfont" << endl;

    outfil << (outer[0].x() + 10) << " mm "
           << (outer[1].y() + 12) << " mm moveto" << endl;

    outfil << "(Template for " << horn_gain << " dBi horn for "
	   <<  int(freq) << " MHz) show" << endl;

    outfil << "-10.0 mm " << (outer[1].y() - 7) << " mm moveto" << endl;

    outfil << "(E-plane) show" << endl;
    outfil << "-10.0 mm " << (outer[1].y() - 14) << " mm moveto" << endl;

    outfil << " (N1BWT 1994) show " << endl;

// and horn outline for single page

    if (pages == 1)
    {
        outfil << "newpath" << endl;
        outfil << PS_moveto() << "\n";
            
        for (int i = 0; i <= 4; i++)
            outfil << PS_lineto(outer[i]) << "\n";

        for ( i = 4; i >= 0; i--)
            outfil << PS_lineto(inner[i]) << "\n";

        outfil << "stroke" << endl;

        outfil << "%" << endl;

// and fold lines 

        outfil << "% fold lines" << endl;


        for ( i = 0; i <= 5; i++)
        {
            outfil << "newpath" << endl;
            outfil << "[10 10] 0 setdash\n" ;
            outfil << PS_moveto(inner[i]) << "\n";
            outfil << PS_lineto(outer[i]) << "\n";
            outfil << "stroke" << endl;
        }

        outfil << "% finish solder flap" << endl;

        outfil << "newpath" << endl;
        outfil << "[10 10] 0 setdash\n" ;
        outfil << PS_moveto(inner[4]) << "\n";
        outfil << PS_lineto(inner[5]) << "\n";
        outfil << "stroke" << endl;

        outfil << "newpath" << endl;
        outfil << "[10 10] 0 setdash\n" ;
        outfil << PS_moveto(outer[4]) << "\n";
        outfil << PS_lineto(outer[5]) << "\n";
        outfil << "stroke" << endl;

        outfil << "% END HORN ***************************************" << endl;
        outfil << "grestore" << endl;

        outfil << "showpage" << endl;

    }
    
    else        // two pages
    {

// E-plane side: coordinates same as single page

outfil << "% E-plane side of HORN ***************************************\n";


        outfil << "newpath" << endl;
        outfil << PS_moveto() << "\n";
            
        for (int i = 0; i <= 1; i++)
            outfil << PS_lineto(outer[i]) << "\n";

        for ( i = 1; i >= 0; i--)
            outfil << PS_lineto(inner[i]) << "\n";

        outfil << "stroke" << endl;

        outfil << "%" << endl;
    
        outfil << "% END E-plane side of HORN ****************************\n";
        outfil << "grestore\n";

        outfil << "showpage\n";
        outfil << "%\n";

// H-plane side: new coordinates

        outfil << "gsave" << endl;
outfil << "% H-plane side of HORN ***************************************\n";

// center large aperture as well as we can


        if (H_aperture > 190)
            offset = Point((0.5 * (215.0 - H_aperture) - outer[0].x()),
                            (254 - outer[0].y()));
                            
        outfil << "%" << endl;
    
    outfil << offset.x() << " mm " << offset.y() << " mm translate" << endl;
    
        outfil << "%" << endl;

        inner[1] = inner[0] + Point(H_wg,0);

        outer[0] = Point(-(0.5 * (H_aperture - H_wg)),H_side);

        outer[1] = outer[0] + Point(H_aperture,0);

        outfil << "newpath" << endl;
        outfil << PS_moveto() << "\n";
            
        for ( i = 0; i <= 1; i++)
            outfil << PS_lineto(outer[i]) << "\n";

        for ( i = 1; i >= 0; i--)
            outfil << PS_lineto(inner[i]) << "\n";

        outfil << "stroke" << endl;

        outfil << "%" << endl;

        outfil << "% END E-plane side of HORN ****************************\n";


// text for H-plane side

        outfil << "/Helvetica-Bold findfont" << endl;
        outfil << "14 scalefont" << endl;
        outfil << "setfont" << endl;
    
        outfil << (outer[0].x() + 10) << " mm "
               << (outer[1].y() + 12) << " mm moveto" << endl;
    
        outfil << "(Template for " << horn_gain << " dBi horn for "
               <<  int(freq) << " MHz) show" << endl;
    
        outfil << "-10.0 mm " << (outer[1].y() - 7) << " mm moveto" << endl;
    
        outfil << "(H-plane) show" << endl;
        outfil << "-10.0 mm " << (outer[1].y() - 14) << " mm moveto" << endl;
    
        outfil << " (N1BWT 1994) show " << endl;

        outfil << "grestore\n";
        outfil << "showpage\n";
        outfil << "%\n";
        
        outfil.close();
    }
}    


//






//* predeclare entry points for phase center routines */
    double Phase_center(char, float, float, float, float);
    double H_Phase_Center(double, double, double, double&, double&);
    double E_Phase_Center(double, double, double, double&, double&);
    double fresnel_sin(double, double&, double&);
    double fresnel_cos(double, double&, double&);


void Antenna::horn_phasecenters()
{

/* get the phase centers and decide what to do about them */

    E_phase_center = Phase_center('E', E_aperture, E_wg, 
                                    wavelength, axial_length);
    H_phase_center = Phase_center('H', H_aperture, H_wg, 
                                    wavelength, axial_length);

//    cout << "E-Plane phase center is "
//	 << setprecision(2) << E_phase_center
//	 << " wavelengths inside horn mouth\n";
//    cout << "H-Plane phase center is "
//         << setprecision(2) << H_phase_center
//         << " wavelengths inside horn mouth\n\n"; 

}

double Pi, Pi2, Pi2sq;     // /* precomputed values for Pi and Pi/2

#define	ERRMARGIN	2.0E-3
#define	MAX(a,b)	(a > b ? a : b)
#define	MIN(a,b)	(a < b ? a : b)


double Phase_center(char plane, float aperture, float guide, 
                    float wavelength, float Horn_length)
//    /* convert to wavelengths for phase center routine */
//    char plane ;        /* E or H plane */
//    double aperture ;   /* aperture dimension in mm */
//    double guide ;      /* waveguide dimension in mm */
    
{
    double side_length ;    /* length of horn side in wavelengths */
    double Ap_lambda ;      /* aperture in wavelengths */
    double v ;              /* sqrt(side_length/2.0) */
    double dmin, dmax ;     /* limits check on phase center */
    double PC ;             /* phase center in wavelengths */

    Ap_lambda = aperture / wavelength ;
    side_length = sqrt(Horn_length*Horn_length +
			 ((aperture - guide)/2.0)*((aperture - guide)/2.0));
    side_length = side_length / wavelength ;
    v = sqrt(side_length/2.0);


    if (plane == 'E')
        PC = E_Phase_Center( Ap_lambda, v, side_length, dmin, dmax);
    else 
        PC = H_Phase_Center( Ap_lambda, v, side_length, dmin, dmax);

/* test for bad phase center */

    if ( fabs(dmax - dmin) != (fabs(dmax) - fabs(dmin)))
        cerr << " \n\nERROR: BUG in "
             << plane << "_Phase_Center function\n\n";// call Matt\n";
    else if ((fabs(dmax - dmin) > (PC /16.0)) && (Ap_lambda > 1.0))
        cerr << " \n\nEXCESSIVE ERROR in "
             << plane << "_Phase_Center function \n\n";
        
    return (PC);
}
//












    
 
/* 

 */

    
///* phase_center.c     Matt Reilly 10/28/91 */
// converted for C++  2/17/94 pcw
///* input data:  a       aperture in wavelengths
//                l       length of horn side in wavelengths
//                v       sqrt(l/2) in wavelengths
//                lv      low de value  - error check
//                hv      high de value -  "     "
//                
//    E_Phase_Center returns distance in wavelengths from mouth of horn to
//        E-plane phase center inside horn.
//
//    H_Phase_Center returns distance in wavelengths from mouth of horn to
//        H-plane phase center inside horn.
//
//*/

void phase_init()
{
    ph_init = 0;
}

double epc(double l, double q, double c2, double s2, double C, double S)
//double l, q, c2, s2, C, S;
{
    return l * ( 1 - q * (c2 * C + s2 * S) / (C * C + S * S));
}

double E_Phase_Center(double a, double v, double l, double &lv, double &hv)
//double a, v, l;
//double * lv, * hv;
{
    double q, c2, s2, fc, fcl, fch, fs, fsl, fsh, res;
    
    if(ph_init) phase_init();
    
    q = a / (2 * v);
    
    fc = fresnel_cos(q, fcl, fch);
    fs = fresnel_sin(q, fsl, fsh);
    
    c2 = cos(Pi2 * q * q);
    s2 = sin(Pi2 * q * q);
    
    res = epc(l, q, c2, s2, fc, fs);
    hv = epc(l, q, c2, s2, fch, fsh);
    lv = epc(l, q, c2, s2, fcl, fsl);
    
    return res;
}


double hpc(double l, double r, double t, double c, double s)
//double l, r, t, c, s;
{
    return l * (1 + (((r * c) + (t * s)) / ((c * c) + (s * s))));
}

double H_Phase_Center(double a, double v, double l, double &lv, double &hv)
//double a, v, l;
//double * lv, * hv;
{
    double U, W, CD, CDH, CDL, SD, SDH, SDL, R, T, u2, w2;
    double su, sul, suh, sw, swl, swh;
    double cu, cul, cuh, cw, cwl, cwh;
    double res;


    if(ph_init) phase_init();
        
    U = v / a + a / (2 * v);
    W = v / a - a / (2 * v);
    u2 = Pi2 * U * U;
    w2 = Pi2 * W * W;
    
    
    R = W * cos(u2) - U * cos(w2);
    T = U * sin(w2) - W * sin(u2);
    
    su = fresnel_sin(U, sul, suh);
    sw = fresnel_sin(W, swl, swh);
    cu = fresnel_cos(U, cul, cuh);
    cw = fresnel_cos(W, cwl, cwh);

    CD = cu - cw;
    CDL = cul - cwh;
    CDH = cuh - cwl;    

    SD = - su + sw;
    SDL = - suh + swl;
    SDH = - sul + swh;    
    
    res = hpc(l, R, T, CD, SD);
    hv = hpc(l, R, T, CDH, SDH);
    lv = hpc(l, R, T, CDL, SDL);

    return res;

}

///* 




///* fresnel_integrals.c         Matt Reilly 10/28/91 */
// converted for C++  2/17/94 pcw
//
///*  approximation to Fresnel sine and Fresnel cosine needed to calculate
//    phase centers of horn 
//    
//*/
//
void fresnel_init()
{
    Pi = atan2(0.0, -1.0);
    Pi2 = Pi / 2.0;
    Pi2sq = Pi2 * Pi2;
    lib_init = 0;
}

double fresnel_f(double z)
//double z;
{
    /* from stegun and abramowitz pp 301 and 302. */
    double sum;
    
    sum = (1.0 + 0.926 * z) / (2.0 + 1.792 * z + 3.104 * z * z);
    
    return sum;
}

double fresnel_g(double z)
//double z;
{
    /* from stegun and abramowitz pp 301 and 302. */
    double sum;
    
    sum = 1.0 / (2.0 + 4.142 * z + 3.492 * z * z + 6.670 * z * z * z);
    
    return sum;
}

double fresnel_cos(double z, double &low, double &high)
//double z, *low, *high;
{
    int n;   /* series index */
    double sum, x, f, g, c, s, fs1, fs2, gc1, gc2, sign;
		 
    sign = 1.0;
    if(lib_init)   fresnel_init();
    
    if (z < 0) 
    {
        z = -z; 
	sign = - 1.0;
    }
    else sign = 1.0;
    
    x = Pi2 * z * z;
    f = fresnel_f(z);
    g = fresnel_g(z);
    c = cos(x);
    s = sin(x);
    

    
    sum = 0.5 + 
	  f * s - 
	  g * c;

    fs1 = (f + ERRMARGIN) * s;
    fs2 = (f - ERRMARGIN) * s;
    gc1 = (g + ERRMARGIN) * c;
    gc2 = (g - ERRMARGIN) * c;

    if(sign > 0)
    {
    	low = (0.5 +  MIN(fs1, fs2) - MAX(gc1, gc2));
    	high = (0.5 + MAX(fs1, fs2) - MIN(gc1, gc2));
    }
    else
    {
        low = -1.0 * (0.5 + MAX(fs1, fs2) - MIN(gc1, gc2));
    	high = -1.0 * (0.5 +  MIN(fs1, fs2) - MAX(gc1, gc2));
    }


    return sum * sign;
}


double fresnel_sin(double z, double &low, double &high)
//double z, *low, *high;
{
    int n;   /* series index */
    double sum, x, err, f, g, c, s, fc1, fc2, gs1, gs2, sign;
		 
		 
    if(lib_init)   fresnel_init();
    

    if (z < 0) 
    {
        z = -z; 
	sign = - 1.0;
    }
    else sign = 1.0;
    
    x = Pi2 * z * z;
    f = fresnel_f(z);
    g = fresnel_g(z);
    c = cos(x);
    s = sin(x);
    
    sum = 0.5 - 
	  f * c - 
	  g * s;

    fc1 = (f + ERRMARGIN) * c;
    fc2 = (f - ERRMARGIN) * c;
    gs1 = (g + ERRMARGIN) * s;
    gs2 = (g - ERRMARGIN) * s;
    
    if(sign > 0.0)
    {
        low = (0.5 - MAX(fc1, fc2) - MAX(gs1, gs2));
    	high = (0.5 - MIN(fc1, fc2) - MIN(gs1, gs2));
    }
    else
    {
        low = -1.0 * (0.5 - MIN(fc1, fc2) - MIN(gs1, gs2));
    	high = -1.0 * (0.5 - MAX(fc1, fc2) - MAX(gs1, gs2));
    }

    return sign * sum;
}
//










//forward declaration
double curve(double , int , double ); 


Lens::Lens() : Antenna()                   //constructor
{
    Lens_diam = 0.0;     
    spacing = 0.0  ;
    width_increase = 0.0;
    N = 0 ;
    BW_E_horn = 0.0 ;
    lens_focal_length = 0.0 ;     
    Gain_lens = 0.0 ;
    index = 0.0 ; 
    compensate = 0;
    step = 0.0;
    lens_to_horn = 0.0;
}    

void Lens::dimension()
//assumes that horn variables are available
{
    int bad = 1;
    while (bad >= 1 )
    {
	cout << "Enter approximate lens diameter in " 
             << units << ": ";
	cout.flush();
	cin >> Lens_diam  ;
        Lens_diam *= metric;
        
	cout << "Enter EXACT lens plate spacing in "
             << units << ": ";
	cout.flush();
	cin >> spacing ;
        spacing *= metric;

	if (spacing < (wavelength/2.0))
        {
	    cout << " ERROR - lens plate spacing below cutoff \n";
            bad++ ;
            if (bad >= 3) break;
        }
        
	else if (spacing >= wavelength)
            {
	      cout << " ERROR - lens plate spacing greater than wavelength \n";
              bad++ ;
              if (bad >= 3) break;
            }
            else bad = 0;
    }
}

void Lens::calculate()
{

//* make lens diameter an even number of spacings */

    N = (int)(Lens_diam /(2.0 * spacing)) ;
    if ( (Lens_diam - (2.0 * N * spacing)) > (spacing * 0.5 ))
        ++N;

    Lens_diam = 2.0 * N * spacing ;
    

    Gain_lens = (Lens_diam / E_aperture) * (Lens_diam / E_aperture) ;
        /* gain over bare horn */    

    Gain_lens =  10.0 * log10( Gain_lens );

    index = sqrt( 1.0 - ( wavelength / (2.0 * spacing)) *
			( wavelength / (2.0 * spacing))) ;
        /* effective index of refraction */

    BW_E_horn = (56.0 / ( E_aperture / wavelength ))/crad ; /* Kraus, p. 380 */
        
    lens_focal_length = Lens_diam / (2.0 * tan(BW_E_horn / 2.0)) ; /* E-plane*/
        
/* NOTE: if the E and H focal planes are different, then we can compensate
         by making the lens have a different focal length in each plane;
         for now, it will be symmetrical using the E-plane focal length.
         The unknown is the distance from the horn to the lens, which is
         the lens focal length minus the horn_focal_point distance. */
        
     
//    cout << " This is an f over "
//         << (lens_focal_length/Lens_diam) << " lens \n";
//    cout << " using " << (2*N+1)
//    	 << " lens plates for an actual diameter of " 
//	 << Lens_diam  << " mm \n\n";
//    cout << "\n";
//    cout << " Estimated gain of lens is " << Gain_lens
//	 << " dB added to horn gain \n\n ":

}

void Lens::phasecheck()
{

/* get the phase centers and decide what to do about them */

//    Horn_E_phase_center = Phase_center('E', Ap_E_horn, guide_E);
//    Horn_H_phase_center = Phase_center('H', Ap_H_horn, guide_H);

    cout << "E-Plane phase center is " << E_phase_center
         << " wavelengths inside horn mouth\n";
    cout << "H-Plane phase center is " << H_phase_center
         << " wavelengths inside horn mouth\n\n";
    

    if (fabs(E_phase_center - H_phase_center)>0.062 /*wavelengths*/)
    {
        compensate = 1 ;
        cout << "     Phase correction may be in order,\n "
             << "     or you might want a different horn\n\n";
    }
    else compensate = 0 ;
}


//void Lens::plates()     // calculate radius of individual plates
//{
//    for (int i = 1 ; i <= N; i++ )
//    {
//        radius[i] = curve(radius[0],i) ;
//        radius_2[i] = curve(radius_2[0],i);
//        width_increase = radius[0] - radius[i] ;
//    }
//}

double curve(double radius, int i, double spacing)
{
    double theta ;       /* angle from axis */
    double radius_N ;    /* radius of plate N */

    theta = asin(( (double)i * spacing) / radius) ;

    radius_N = radius * cos(theta) ;

    return (radius_N) ;
}


void Lens::write_data(ostream& ostrm, int interactive)
{
    int i ;
    double FGHz = c / (wavelength * 1000000);

/* now a nice header */

    ostrm << "      Metal plate lens antenna for microwaves  \n\n";
    ostrm << "        ALL dimensions are in millimeters! \n\n" ;

    ostrm << "           ____________     \n";
    ostrm << "          |     E      |    \n";
    ostrm << "          |     ^ -->H |   waveguide  \n";
    ostrm << "          |     |      |    \n";
    ostrm << "           ------------     \n\n";

    ostrm << "  At a center frequency of "
	  << setprecision(3) << FGHz << " GHz\n\n";

    ostrm << "  For a lens with a diameter of "
	  << setprecision(2) << Lens_diam << " mm."
	  << " and a plate spacing of "
	  << setprecision(3) << spacing << " mm. \n\n";

    ostrm << "  Fed by a horn of axial length = "
	  << setprecision(2) << axial_length << " mm, \n"
	  << "      H-plane aperture = " << H_aperture << " mm \n"
	  << "      E-plane aperture = " << E_aperture << " mm \n "
	  << "      and a Gain of " << horn_gain
	  << " dB over isotropic \n\n";

    ostrm << "   E-Plane phase center is " << setprecision(2)
	  << E_phase_center << " wavelengths inside horn mouth\n";
    ostrm << "   H-Plane phase center is " << setprecision(2)
	  << H_phase_center << " wavelengths inside horn mouth\n";

    lens_to_horn = lens_focal_length - (E_phase_center * wavelength);

    ostrm << "  Calculations for an f/" << setprecision(2)
	  << (lens_focal_length/Lens_diam) ;
    ostrm << " lens with a focal length of " << lens_focal_length << " mm. \n"
	  << "      providing an estimated gain of " << Gain_lens
	  << " db over the horn\n\n"
	  << "  Distance from horn mouth to center of lens curve is "
	  << lens_to_horn << " mm.\n";

//    ostrm << "\n\n";
      if(interactive)
      {
	    char pause;
	    cout << "\n\nEnter D to display lens design dimensions: ";
	    cout.flush();
	    cin >> pause;
       }
/* design data */

    double radius[25];
    double radius_2[25];
    int SavedN;
/*    radius[0] = (index - 1.0) * lens_focal_length ; */ /* correct form */
    radius[0] = (1.0 - index ) * lens_focal_length ;      /* positive number */

    radius_2[0] = 2.0 * radius[0] ;

    SavedN = N;
    if (N > 24)
    {
	ostrm << "This is a pretty big lens, so only 50 plates calculated\n";
	N = 24;
    }

    for (i = 1 ; i <= N; i++ )
    {
	radius[i] = curve(radius[0],i,spacing) ;
	radius_2[i] = curve(radius_2[0],i,spacing);
    }

    if (!compensate)
    {

    ostrm << "\n  Radius of Curvature of lens plates starting from "
	  << "center plate\n\n" ;
    ostrm << "   Plate  Single radius   double radius   plate width \n";
    ostrm << "   -----    -------         ---------      -------- \n\n";
    ostrm << setprecision(2)
	  << "   ctr     " << radius[0] << " mm.     " << radius_2[0]
	  << " mm.    min \n";

	for (i = 1 ; i <= N; i++ )
	{
	    width_increase = radius[0] - radius[i] ;

	    ostrm << setprecision(2)
	    << "    " << i << "      " << radius[i] << " mm.     " << radius_2[i]
	    << " mm.    min + " << width_increase << "\n";
	}
    }

    else            /* add column for compensation */
    {

    double c_focal_length = lens_focal_length +
		     (H_phase_center - E_phase_center) * wavelength;

    double c_radius[25];
    c_radius[0] = (1.0 - index) * c_focal_length ;

    ostrm << "\n  Radius of Curvature of lens plates starting from "
	  << "center plate\n" ;
    ostrm << "  with H-Plane phase compensation by shifting plate "
	  << "centers\n";
    ostrm << "  of double-curved lens; positive shift is toward horn.\n\n";

    ostrm << "   Plate  Single radius   Double Radius   Plate Width"
	  << "      Shift\n";
    ostrm << "   -----    -------         ---------     -----------"
	  << "      ----- \n\n";

    ostrm << setprecision(2)
	  << "   ctr     " << radius[0] << " mm.     " << radius_2[0]
	  << " mm.       min \n";

    double shift [25];

	for (i = 1 ; i <= N; i++ )
	{
	    width_increase = radius[0] - radius[i] ;
	    c_radius[i] = curve(c_radius[0],i,spacing);
	    shift[i] = radius[0] - radius[i] - (c_radius[0] - c_radius[i]);

    ostrm << setprecision(2)
	  << "    " << i << "      " << radius[i] << " mm.     " << radius_2[i]
	  << " mm.       min + " << width_increase << "   "
	  << (shift[i]/2.0) << "\n";
	}
    }

    step = wavelength / ( 1.0 - index );

    if ( step < width_increase)
    ostrm << setprecision(2) << "\nFor a Zoned lens plate, "
	  << "step distance = " << step <<" mm.\n\n";

    N = SavedN;
}

void Lens::crude_graphics(ostream& ostrm)
{
/* and a picture to help */


ostrm << "\n\n\n       ****** UNBELIEVABLY CRUDE GRAPHICS ******\n\n";
                                                                                
ostrm << "               Single Curve        Double Curve \n\n";   
                                                                             
ostrm << "                       -----------         -------------\n";
ostrm << "                        ----------          ----------- \n";
ostrm << "        __               ---------    ^      --------- \n";
ostrm << "       |  |               --------    |       -------  \n";
ostrm << "       |  |    <          --------   H|       ---+---  \n";
ostrm << "       |  |    horn       --------    |       -------  \n";
ostrm << "        --               ---------           --------- \n";
ostrm << "       waveguide        ----------          -----------\n";
ostrm << "                       -----------         -------------\n";
ostrm << "                                                 +\n";
ostrm << "                                    centerline   + \n";
ostrm << "                                                 +\n";
ostrm << "                       -----------         -------------\n";
ostrm << "                        \\  LENS  |         \\   LENS    /\n";
ostrm << "       waveguide         \\ PLATE |    ^     \\  PLATE  /\n";
ostrm << "        _____             \\      |    |      \\       / \n";
ostrm << "       |     | <           )<--->|   E|       )<-+->(  \n";
ostrm << "        -----  horn       / min. |    |      /  min. \\ \n";
ostrm << "                      -->/ thick.|       -->/  thick. \\ \n";
ostrm << "                     /  /        |      /  /           \\ \n";
ostrm << "                    /  -----------     /   -------------\n";
ostrm << "                   /                  /                 \n";
ostrm << "                   Radius of curvature            \n\n\n";
}
    





// ****** end of ANTENNA.CPP file



// start of HDL_ANT.CPP 





//forward declarations
int KeyEvent(void);
void menu();
void pdish();
void mlens();
void range();
void newhorn();
void oldhorn();
void nf_corr();
void info();
void postinfo();
void english_units();

int main()
{

    Antenna start;
    start.banner(cout);

    cout << "\n\n\n HDL_ANT - designs, makes calculations, and draws \n"
	 << "           construction templates for microwave \n"
         << "           horn and parabolic dish antennas . \n\n"
	 << " N1BWT 1994 \n\n";
//pause
    char pause;
    cout << "Hit <ENTER> key to continue";
    cout.flush();
    cin.setf(0,ios::skipws);
    cin >> pause;
    cin.setf(1,ios::skipws);

    int k = 0;
    menu();
    
    while ( (k = KeyEvent() ) != CTRL_C)
   {
	if (k)
	{

	   // Mask out the key's scan code.
	   k &= 0x00ff;
	   char mode = (char)k;

	   mode = toupper(mode);
        
	    switch(mode)
            {
                case 'H':   newhorn(); break;    
                case 'E':   oldhorn(); break;
		case 'D':   pdish(); break;
                case 'L':   mlens(); break;
		case 'R':   range(); break;
		case 'C':   nf_corr(); break;
		case 'I':   info(); break;
		case 'P':   postinfo(); break;
                case 'U':   english_units(); break;
                case 'Q':   return 0;
                default: cout  << "Huh?"; break;
	    }

            char pause;
            cout << "\n\nEnter C to continue: ";
	    cout.flush();
	    cin >> pause;

	    menu();
        }             
    }
   return 0;
}
                 
void menu()             // just prints a menu to select from
{                

// clearscreen and reset xy
//        textmode(MONO);
	clrscr();
        gotoxy(1,1);


	cout << "\n"
	  << "Enter first letter for selection \n\n"
          << "   Horn antenna design and template\n"
          << "   Existing horn antenna calculations\n"
          << "   Dish antenna calculations and template\n"
          << "   Lens antenna design\n"
          << "   Range design for antenna measurement\n"
          << "   Corrections for antenna measurements\n"
          << "   Information about HDL_ANT program\n"
          << "   PostScript printing information\n"
          << "   Units: Metric [default] or English\n"
          << "   Quit\n";
}

int KeyEvent(void)
{
    // Check for key press.
    unsigned key = bioskey(1);
    
    // Get key if one is  available.
    if (key) key = bioskey(0);
    
    return key;
    
}


void pdish()
{
    cout << "\n\nPARABOLIC DISH ANTENNA CALCULATIONS AND TEMPLATE GENERATION\n\n";
    Dish d;
    d.get_freq();
    d.d_init();
    d.illumination();
    d.calc_gain();   
    char mode = 'n';
    cout << "\nDo you want to make a PostScript template [Yes or No]?  ";
    cout.flush();
    cin >> mode ;
    mode = toupper(mode);

    if (mode == 'Y') d.PS_template();
}


void newhorn()
{
    cout << "\n\nDESIGN A HORN ANTENNA AND GENERATE A CONSTRUCTION TEMPLATE\n\n";
    char mode = 'Y';
    Horn a;
    a.get_freq();
    a.get_wg();
    if (a.bad_input() >= 3) goto end;
    a.lookup();

    while (mode == 'Y')
    {
        a.simple_gain();
        a.horn_angles();
        a.horn_calcgain();
        a.horn_phasecenters();
	a.write_data(cout,1);

        cout << "\nDo you want to change horn dimensions [Yes or No]?  ";
	cout.flush();
	cin >> mode ;
        mode = toupper(mode);

        if (mode == 'Y') a.horn_dimensions();
    }

    mode = 'n';
    cout << "\nDo you want to make a PostScript template [Yes or No]?  ";
    cout.flush();
    cin >> mode ;
    mode = toupper(mode);

    if (mode == 'Y')
    {
        a.coords();
        a.PS_template();
    }

    end: ;
}


void oldhorn()
{
    cout << "\n\nHORN ANTENNA CALCULATIONS AND CONSTRUCTION TEMPLATE "
         << "GENERATION\n\n";
    char mode = 'n';
    Horn a;
    a.get_freq();
    a.get_wg();
    if (a.bad_input() >= 3) goto end;
    a.horn_dimensions();
    a.simple_gain();
    a.horn_angles();
    a.horn_calcgain();
    a.horn_phasecenters();
    a.write_data(cout,1);

    cout << "\nDo you want to make a PostScript template [Yes or No]?  ";
    cout.flush();
    cin >> mode ;
    mode = toupper(mode);

    if (mode == 'Y')
    {
        a.coords();
        a.PS_template();
    }

    end: ;
}


void mlens()
{
    cout << "\n\nMETAL PLATE LENS ANTENNA DESIGN\n\n";
    char mode = 'n';
    Lens m;
    m.get_freq();
    m.dimension();
    if (m.bad_input() >= 3) goto end;
    m.get_wg();
    m.horn_dimensions();
    m.horn_angles();
    m.horn_calcgain();
    m.calculate();
    m.horn_phasecenters();
    m.phasecheck();
    m.write_data(cout,1);

// set up output file if desired;

    cout << " Save this lens design to file (yes or NO) ";
    cout.flush();
    cin >> mode ;
    mode = toupper(mode);

    if (mode == 'Y')
    {
	char  outname[LINESIZE];
        cout << "\nEnter output filename: " ;
	cout.flush();
	cin >> outname ;
    
        ofstream ostrm(outname);

        m.banner(ostrm);
	m.write_data(ostrm,0);
        m.crude_graphics(ostrm);
    }
    end: ;
}

void range()
{
    cout << "\n\nANTENNA RANGE DESIGN CALCULATIONS\n\n";
    Antenna range;
    range.get_freq();
    range.range_calc();
}

void nf_corr()
{
    cout << "\n\n Correction of PANFI readings to Antenna Gain\n\n";
    Antenna c;
    c.get_freq();
    c.PANFI_corr();
}

void info()
{
    Antenna info;
    info.banner(cout);    
    info.info_listing(cout);

// set up output file if desired;

    char mode;
    cout << " \n\nWrite this information to file (Yes or No) ";
    cout.flush();
    cin >> mode ;
    mode = toupper(mode);

    if (mode == 'Y')
    {
	char outname[LINESIZE];
        cout << "\nEnter output filename: " ;
	cout.flush();
	cin >> outname ;
    
        ofstream ostrm(outname);

        info.banner(ostrm);
        info.info_listing(ostrm);
    }
}

void postinfo()
{
    Antenna pp;
    pp.post_print(cout);

// set up output file if desired;

    char mode;
    cout << " \n\nWrite this information to file (Yes or No) ";
    cout.flush();
    cin >> mode ;
    mode = toupper(mode);

    if (mode == 'Y')
    {
	char outname[LINESIZE];
        cout << "\nEnter output filename: " ;
	cout.flush();
	cin >> outname ;
    
        ofstream ostrm(outname);

        pp.post_print(ostrm);
    }
}

void english_units()
{
    char mode;

    cout << " \n\nCurrently using " << unit_type 
         << " units for dimensions\n";
    
    cout << " \n\nChange to alternate units(Yes or No) ";
    cout.flush();
    cin >> mode ;
    mode = toupper(mode);

    if (mode == 'Y')
    {
        if (metric == 1.0)  // change to English
        {
            metric = 25.4;
	    units = "inches";
	    unit_type = "English";
	}
	else                // change to metric
	{
	    metric = 1.0;
	    units = " mm.  ";
            unit_type = "Metric ";
        }
    }
}
         
    

// end of HDL_ANT.CPP