/*Interdigital Bandpass Filter design program  V 1.01  [Linux version]

Conversion from Basic to C++ by Dale Heatherington, WA4DSY,  5-16-1996
--------------------------------------------------------------------

May 19 2004
Fixed a little error in cheb() pointed out by Federico Massaioli 
and changed some of the includes so this will compile on current
versions of GNU C++ .

-----------------------------------------------------------------------

This version is designed as a CGI program to allow
it to run from a HTTP server where it can be accessed by
Web Browsers. It outputs HTML formated text. Input
is from the Browsers FORM.

The original BASIC  program was published in the January 1985 issue of
Ham Radio magazine on page 12.  The authors are Jerry Hinshaw, N6JH
and Shahrokh Monemzadeh.


Converted from OS/2 to Linux 8-13-97  by Dale Heatherington

Compile with this command:

$g++ -o idbpf idbpf.cpp

Note: The input to this program comes from environment variable "QUERRY_STRING"
which is supplied by the www server.

*/


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

#include <iostream>
#include <iomanip>
#include <fstream>

#include <math.h>
using namespace std;


//These arrays are copied directly from the BASIC code but start at zero instead of 1.

double G[200],C[200], RK[200], AK[200], FR[40],ALOSS[40], A[200], B[200] ;

//--------------
// These globals supply the program input data

int elements;            //number of elements
double ripple;              //ripple if cheby
double FZGC,BWMC,R;      //f center in GHZ, BW in mhz, load resistance in ohms
double H,D,E;            // Height above gnd, rod dia, end spacing
int NFR;                     //number of frequency rejection points in graph
double STP;                  //frequency step size

//--------------

const double PI = 3.14159265;

const char instructions[] =
"Fill out the form below to design a bandpass filter as"
" described in the January 1985 issue of Ham Radio Magazine."
" These filters are generally limited to frequencies above 200 MHZ "
" because their size is slightly longer than 1/4 wavelength. "
" <A href=\"/idbpfdraw.html\">Click here for a drawing of "
" a generic filter.</A>"
" You may also download the "
" <A href=\"/idbpf.txt\"> C++ source</A> for this program";



//--------------------------------------------------------------

double FuncFNRJ(double ta, double b, double c, double d)
{
  return ((b*c)-(ta*d)) / ((c*c) + (d*d));
}

//--------------------------------------------------------------

void cheb(int ele, double rip)
{
  double beta,gamma,c;
  int k;

  c = (2 * rip) / 17.37;
  beta = log((exp(c)+1) / (exp(c) -1));
  gamma = 0.5 * (exp(beta/(2*ele)) - exp(-beta / (2*ele)));
  for(k=0;k<ele;k++)
     {
        A[k] = sin(0.5 * (2 * (k+1) - 1) * PI / ele);
        B[k] = (gamma * gamma) + pow(sin((k+1)*PI/ele) , 2);
     }
  G[0] = 2 * A[0] / gamma;
  for(k=1;k<ele;k++)
     {
        G[k] = 4 * A[k-1] * A[k] / (B[k-1] * G[k-1]);
     }
  int nn = ele / 2;
  int nnn = (ele+1) / 2;
  int s = nnn - nn;
  if ( s <= 0)
     { G[ele] = pow(((exp(beta/2)+1) / (exp(beta/2)-1)) , 2);
        return;
     }
  G[ele] = 1;
  return;

}

//--------------------------------------------------------
void butter(int ele)
{
  const double pov2 = 1.57079633;
  int k;

  for(k=0;k<ele;k++)
     {
        G[k] = 2 * sin(pov2 * (2*(k+1) - 1) / ele);

     }

  G[ele] =  1;
  return;
}

//-----------------------------------------------------------
void Compute(void)
{
  int k,ip,conter=0;
  int start = -NFR / 2;
  int end = NFR / 2;
  double bw3gc,bwrgc, bw3,b,ca;
  int ele = elements;

//line 190
  for(ip=start;ip<end;ip++)
     {
        FR[conter++] = FZGC + (STP * 0.001 * ip);
     }

  double F1 = FZGC - 0.0005 * BWMC;
  double F2 = FZGC + 0.0005 * BWMC;

  if (ripple == 0.0)
     {
        bw3gc = F2 - F1;
        bwrgc = 0;
        bw3 = 1;
        butter(ele);
     }else
        {
           b = 1 / sqrt(pow(10,0.1*ripple)-1);
           ca = log(b + sqrt(b * b - 1)) / ele;
           bw3 = (exp(ca) + exp(-ca)) / 2;
           cheb(ele,ripple);
           bwrgc = F2 - F1;
           bw3gc = bwrgc * bw3;
         }

     double w = 2 * (F2 - F1) / (F2 + F1);
     double qf = FZGC / bw3gc;
     int nfm = ele -1;
     double qwvl = 11.8028 / (4 * FZGC);
//line 440
     for(k = 0;k<nfm;k++)
        {
           AK[k] = 1 / (bw3 * sqrt(G[k] * G[k+1]));
           RK[k] = AK[k] / qf;
        }

     double ako = G[0] * bw3;
     AK[ele-1] = ako;
     AK[ele] = 0;
     double qs = G[0] * bw3 * qf;
     double canh = (exp(2*PI*E/H)-1) / (exp(2*PI*E/H)+1);
     double zm = 59.9585 * log(4*H / (PI*D));
     double ze = 59.9585 * log(canh * H * 4 / (PI * D));
     double rkm = RK[0] * sqrt(zm/ze);
     double z = PI * D / (2 * H);
     double coth = (exp(z) + 1) / (exp(z) -1);
     double y = PI * rkm / 4;
     double t = pow(coth,y);
//line 600
     C[0] = (H / PI) * log((t+1) / (t-1));
     int mfl = ele -2;
     for(k=2;k<=mfl;k++)
        {
           double y = PI * RK[k-1] / 4;
           double t = pow(coth,y);
           C[k-1] = (H / PI) * log((t+1) / (t-1));
        }
     C[ele-2] = C[0];
     double x = sqrt(PI * R / (4 * ze * qs));
     double aq =  2 * qwvl * atan(x / sqrt(1 - x * x)) / PI;
     double qu = 2200 * H * sqrt(FZGC);
     double sumg = 0;
//line 740
     for(k=0;k<ele;k++) sumg += G[k];
     double bloss = 4.34 * FZGC * sumg / (qu * (F2 - F1));
     double delay = sumg / (2 * PI * (F2 - F1));

     cout.setf(ios::showpoint | ios::fixed);

     cout  << setprecision(3)
           << "Design data for "
           << ele
           << " section interdigital bandpass filter.\n"
           << "Bandpass ripple = "
           << ripple
           << " DB\n"
           << "Center frequency = "
           << FZGC * 1000
           << "MHZ\n"
           << "Cutoff frequency = "
           << F1*1000 << " and " << F2*1000 << " MHZ\n"
           << "Ripple BW = "
           << bwrgc * 1000 << " MHZ\n"
           << "3 DB BW = "
           << bw3gc*1000 << " MHZ\n"
           << "Fractional BW = " << w*1000 << " MHZ\n"
           << "Filter Q = " << qf << endl
           << "Est QU = " << qu << endl
           << "Loss Based on this QU = " << bloss << " DB\n"
           << "Delay at band center = " << delay << " NS\n" ;


     cout  << setw(71) << "Freq"
           << setw(9) << "DB loss" << endl;

	  int i;
     for(i=0;i<NFR;i++)
        {
           double nfn = (2 * (FR[i] - FZGC)) / (w * FZGC);
           if (nfn < 0.0) nfn = -nfn;  //abs value of nfn

           if (ripple == 0.0)
              {
                 ALOSS[i] = 10 * log((double)1 + pow(nfn,2*ele)) / log((double)10);
              }
              else
                 {
                    if (nfn < 1.0) nfn = 1.0;
                    double ang = ele * log(nfn + sqrt(nfn * nfn - 1));
                    double yak = 0.5 * (exp(ang) + exp(-ang));
                    ALOSS[i] = 10 *  log((double)1 + (pow(10,(0.1*ripple)) - 1) * yak * yak) / log((double)10);
                 }
            int Aloss ;
            if (ALOSS[i] > 65) Aloss = 65 ; else Aloss = (int)ALOSS[i];
            double fr = FR[i] * 1000;
            double Alos = ALOSS[i];


// Print loss vs frequency graph.

            cout << setw(Aloss+1) << "*"
                 << setprecision(1)
                 << setw(72-Aloss) << fr
                 << setw(7) << Alos << endl;
        }

//line 1060

     double wo = 2 * PI * FZGC * 1e+9;
     double f = D / H;
     double cf = (-0.0000422 + 0.0857397 * f + 0.0067853 * f * f - 9.092165e-2
                    * pow(f,3) + 0.169088 * pow(f,4)) * PI * H * 2.54;

     double ww   = wo * 1e-12;
     double b2   = PI * aq / ( 2 * qwvl);
     double gg   = 1 / R;
     double bb   = -cos(b2) / (ze * sin(b2));
     double el1  = 0.8 * qwvl;

     double ang  = el1 * PI / (2 * qwvl);
     double b1   = ang - b2;
     double yl   = -cos(ang) / (zm * sin(ang));
     double cp   = ww * (cf + 0.17655 * D * D / (qwvl - el1));
     double y1   = cp + yl;
     double el2  = 0.87 * qwvl;

            ang  = el2 * PI / ( 2* qwvl);
     double b4   = ang - b2;
            yl   = -cos(ang) / (zm * sin(ang));
     double cd   = ww * (cf + 0.17655 * D * D / (qwvl - el2));
     double y2   = cd + yl;
     double el3  = 0.95 * qwvl;

            ang  = el3 * PI / (2 * qwvl);
     double  b5  = ang - b2;
             yl  = -cos(ang) / ( zm * sin(ang));
     double  cq  = ww * (cf + 0.17655 * D * D / (qwvl - el3));
     double  y3  = cq + yl;
     double elem = y3 * y2 * el1 / ((y1 - y2) * (y1-y3)) + y1 * y3 * el2
                    / ((y2-y1)*(y2-y3)) + y1 * y2 * el3 / ((y3-y1) * (y3-y2));

     double tann = sin(b1) / cos(b1);
             yl  = FuncFNRJ(gg, bb+tann/ze, 1-ze*bb*tann, ze*gg*tann);
             y1  = cp + yl;
            tann = sin(b4) / cos(b4);
              yl = FuncFNRJ(gg, bb+tann/ze, 1-ze*bb*tann, ze*gg*tann);
             y2  = cd + yl;
            tann = sin(b5) / cos(b5);
              yl = FuncFNRJ(gg, bb+tann/ze, 1-ze*bb*tann, ze*gg*tann);
             y3  = cq + yl;

     double eleq = y3 * y2 * el1 / ((y1 - y2) * (y1 - y3)) + y1 * y3 * el2
                   / ((y2 - y1) * (y2 - y3)) + y1 * y2 * el3 / ((y3 - y1) * (y3 - y2));

     cout  << setprecision(3)
           << setw(30) << "Quarter wavelength = " << qwvl << " inches\n"
           << setw(30) << "Length of interior elements = " << elem << " inches\n"
           << setw(30) << "Length of end elements = " << eleq << " inches\n"
           << setw(30) << "Ground plane space = " << H << " inches\n"
           << setw(30) << "Rod diameter = " << D << " inches\n\n"
           <<  "End plates are " << E << " inches from center of end rod.\n"
           << "Tap external lines up " << aq << " inches from shorted end.\n"
           << "Line impedances: End rod: " << ze << ", other: " << zm
           <<  " Ext. lines: " << R << " Ohms\n";


int Tab = 10;

    cout   << "Dimensions (inches)\n"
           << setw(Tab) << "El. No."
           << setw(Tab) << "End to C"
           << setw(Tab) << "C to C"
           << setw(Tab) << "G[k]"
           << setw(Tab) << "Q/Coup\n" ;

    double dom = E;
    int goo = 1;

    cout   << setw(Tab) << "0" << endl
           << setw(Tab) << "1"
           << setw(Tab) << E
           << setw(Tab) << C[0]
           << setw(Tab) << G[0]
           << setw(Tab) << AK[0]
           << endl;


  for(i=0;i<nfm;i++)
       {
           int L = i + 1;
           dom = dom + C[i];
           cout  << setprecision(3)
                 << setw(Tab) << L+1 << setw(Tab) << dom
                 << setw(Tab) << C[L]
                 << setw(Tab) << G[L]
                 << setw(Tab) << AK[L]
                 << endl;
       }

    cout   << setw(Tab) << ele+1
           << setw(Tab) << dom+E
           << endl;

    cout   << endl
           << "Box inside dimensions: " << qwvl
           << "H x " << dom+E
           << "W x " << H
           << "T  inches\n";


}
//------------------------------------------------------------
char *Validate(void)
{

  if ((elements < 2) || (elements > 12)) return "Elements";
  if ((ripple < 0) || (ripple > 4)) return "Ripple";
  if (FZGC <= 0) return "Frequency";
  if (BWMC <= 0) return "Bandwidth";
  if (R < 1)    return "Impedance";
  if ((H <= 0) || (H <= D)) return "Ground plane space";
  if (D <= 0)  return "Rod Dia";
  if (E <= (D/2)) return "End plate to rod";
  if (NFR < 3) return "Points to graph";
  if (STP <= 0) return "Step size";

  return NULL;

}
//------------------------------------------------------------
//query_string format :  name=value&name2=value2&name3=value3&.....name10=value10

void ProcQuery(char *cp)
{
  char* arg;
  char* name;
  float val=0.0;

	char *s = strtok(cp,"&");       //s = name,value pair into s

	while(s != NULL)
    {
	  	  name = s;
		  arg = strchr(s,'=');
		  *arg = (char)NULL;
		  arg++;

		  sscanf(arg,"%f",&val);

        if (strcmp(name,"ele") == 0)   elements = (int)val;
        if (strcmp(name,"rip") == 0)   ripple = (double)val;
        if (strcmp(name,"fzmc") == 0)  FZGC = (double)val/1000;
        if (strcmp(name,"bwmc") == 0)  BWMC = (double)val;
        if (strcmp(name,"R") == 0)     R = (double)val;
        if (strcmp(name,"H") == 0)     H = (double)val;
        if (strcmp(name,"D") == 0)     D = val;
        if (strcmp(name,"E" ) == 0)   E =(double)val;
        if (strcmp(name,"NFR") == 0)   NFR = (int)val;
        if (strcmp(name,"STP") == 0)   STP = (double)val;

			s = strtok(NULL,"&");		 //get next name,value pair
    }


}

//-------------------------------------------------------------
void SendForm(void)
{
     cout  << "<FORM ACTION=\"idbpf\" METHOD=GET>"
        << "<p><h2>Design a custom Interdigital Bandpass Filter</h2>"
        << instructions
        << "<hr><br><INPUT TYPE=\"text\" NAME=\"ele\" SIZE=5 > Number of elements"
        << "<br> <INPUT TYPE=\"text\" NAME=\"rip\" SIZE=5 value=\"0\"> Ripple in DB "
         << "<br> <INPUT TYPE=\"text\" NAME=\"fzmc\" SIZE=5 > Center Frequency in MHZ "
        << "<br> <INPUT TYPE=\"text\" NAME=\"bwmc\" SIZE=5 > Bandwidth in MHZ"
         << "<br> <INPUT TYPE=\"text\" NAME=\"R\" SIZE=5 value=\"50\"> I/O Impedance in OHMS"
         << "<br> <INPUT TYPE=\"text\" NAME=\"H\" SIZE=5 > Ground plane space in INCHES"
        << "<br> <INPUT TYPE=\"text\" NAME=\"D\" SIZE=5 > Rod diameter in INCHES"
         << "<br> <INPUT TYPE=\"text\" NAME=\"E\" SIZE=5 > End plate to rod in INCHES"
         << "<br> <INPUT TYPE=\"text\" NAME=\"NFR\" SIZE=5 value=\"30\"> Points to graph (max = 40) "
         << "<br> <INPUT TYPE=\"text\" NAME=\"STP\" SIZE=5 value=\"2\"> Frequency step size in MHZ"
        << "<p> <INPUT TYPE=\"submit\">  <INPUT TYPE=\"reset\">"
        << "<hr>"
        << "</FORM></BODY></HTML>"
        << endl;

}



//-----------------------------------------------------------

//This is a http CGI program.  Input is from the environment variable QUERY_STRING.

int main(int argc, char *argv[])
{

  char *cp;

  char *query_string = getenv("QUERY_STRING");
  if (query_string == NULL) query_string = "";

	cout  << "Content-type: text/html"
        << endl
        << endl
        << "<HTML><HEAD><TITLE>Interdigital Bandpass Filter Designer</TITLE></HEAD>"
        << "<BODY BGCOLOR=\"#FFFFFF\" >" ;



  if (strlen(query_string) == 0)
     {
        SendForm();

     }
  else
     {
        ProcQuery(query_string);       //Parse the query string
         cp = Validate();              //validate user input
        if ( cp == NULL)
           {
              cout << "<pre>\n";
              Compute();               //Compute filter dimensions
              cout << "</pre>\n";
           }
           else
              cout  << "<br><H3>Invalid data entered: " << cp << " error</H3>\n";
     }


   cout << "</body></html>" << flush;
  return 0;
}