program fildisp;

{               Filter Response & Stores Data to Disk        }

type     XYData = record
                        XX,YY: real;
                  end;
         vector = array[1..151] of real;
         vector1= array[0..20] of real;
         vector20 = array[0..20] of real;
         Plotarray = array[1..100,1..2] of real;

var      A              : Plotarray;
         graph_title    : string[15];
         N              : integer;
         Fr,mag         : vector;
         wp_lo,wp_hi    : real;
         wo_geom        : real;
         w_beg,w_end    : real;
         order          : integer;
         eta            : real;
         w_inc          : real;
         km             : integer;
         w              : real;
         i,j,k,m        : integer;
         wp             : real;
         p,p_2          : real;
         B,r            : real;
         ak,bk,g        : vector1;
         gk_sum         : real;
         L_dB           : real;
         Qu             : real;
         T_do           : real;
         Td_sum         : real;
         Peak           : real;
         Loss,Loss_wp   : real;
         phi            : real;
         Anz,Ans1,Ans2  : char;
         T_Rhs,T_Lhs    : real;
         Value          : string[20];
         Const_2,Const_1: real;
         Pole_Re        : vector20;
         Pole_Im        : vector20;
         Time           : vector;
         t1,t2,t3       : real;
         Odd_Order      : Boolean;
         Ang,Ang_1      : real;
         xmax,xmin,ymax,ymin: real;
         Title              : string[15];
         file_name          : file;

label 500;
{                              Math Library Functions                       }
{                                   Version  2.0                            }
{                                                                           }
Function Tan( X: Real): real;
begin
    Tan:= Sin(X)/Cos(X);
end;
{ ///////////////////////////}

Function Sign ( X:Real ): real;
    var y: real;
    begin
         if X < 0 then Sign:= -1.0 else Sign:= 1.0;
    end;
{ ///////////////////////////}

Function Atan4( Y,X:Real ): Real ;  {  4 Quadrant ArcTan  }
    const Pi = 3.141592654;
    var T1 : real;
    begin
         if abs(X) > 1.0E-8 * ABS(Y) then
              begin
                   T1:= ArcTan( Y/X );
                   if X < 0 then T1:= T1 + Pi;
              end
              else T1:= 0.5*Pi*Sign ( Y );
         Atan4:= T1;
    end;
{ ///////////////////////////}

Function Log10( X: real): real;
    begin
         Log10:= 0.434294482*Ln( X );
    end;
{ /////////////////////////// }

Function Cosh ( X: real ): real;
    var temp: real;
    begin
         temp:= exp( X );
         Cosh:= 0.5* (temp + 1/temp );
    end;
{ /////////////////////////// }

Function Sinh ( X:real ): real;
    var temp : real;
    begin
         temp:= exp( X );
         Sinh:= 0.5* (temp - 1/temp);
    end;
{ /////////////////////////// }

Function ArcSinh( X: real) : real;
    begin
         ArcSinh:= Ln ( X + sqrt( sqr(X) + 1 ));
    end;
{ /////////////////////////// }

Function XtoY ( X,Y: real): real;
    begin
         XtoY:= exp( Y*Ln(X));
    end;
{ /////////////////////////// }

Function Tanh( X: real): real;
    var temp,temp2: real;
    begin
         temp:= exp( X );
         temp2:= 1/temp;
         Tanh:= ( temp -temp2 )/( temp + temp2 );
    end;
{ /////////////////////////// }

Function ArcCosh( X: real): real;
    begin
         ArcCosh:= Ln ( X + sqrt( sqr(X) - 1));
    end;
{ /////////////////////////// }

Function ArcCos( X: real): real;
    begin
         ArcCos:= arctan( sqrt( 1 - sqr(X))/X );
    end;
{ /////////////////////////// }

Procedure Write_to_File( X,Y: vector);

var j              : integer;
    write_file     : file of XYData;
    read_file      : file of XYData;
    file_rec       : XYData;
    file_name      : string[10];
    file_size      : integer;
    OK             : Boolean;
    xyz            : boolean;

label 1;

begin
clrscr;
gotoxy(10,10);
write('Store Inputed Data to Disk  y / n  ?  ');
readln(Ans1);
if Ans1 = 'Y' then
    begin
    clrscr;
    OK:= False;
    gotoxy(10,10);
    write('Write Data to What File  ?  ');
    readln(file_name);
    Assign(write_file,file_name);
    {$I-} Reset( write_file)  {$I+};
    OK:= (IOResult = 0);   { Determine if Write-to File Already Exists }
    if OK then
         begin
         clrscr;
         gotoxy(10,10);
         write('File Already Exists . . .  Overwrite  y / n  ? ');
         readln(Ans2);
         if Ans2 <> 'y' then goto 1;
         end;
    Assign(write_file,file_name);
    Rewrite(write_file);
    with file_rec do
         begin
         for j:= 1 to km do   { Writes Data to Disk File }
              begin
              XX:= X[j];
              YY:= Y[j];
              write(write_file,file_rec);
              end;
         end;
    close(write_file);
 end;
 1: delay(10);
 end;


begin
500: clrscr;
gotoxy(1,5);
write('     Enter    Low,    Hi    Ripple Freqs,   MHz  ');
readln(wp_lo,wp_hi);
write('     Enter    Start   Stop   Freqs    MHz  ');
readln(w_beg,w_end);
write('     Enter   Freq   Increment    MHz  ');
readln(w_inc);
write('     Enter   dB   Passband   Ripple  ');
readln(r);
write('     Enter   Order  of  Filter;    Unloaded Q   ');
readln(order,Qu);
eta:= sqrt(XtoY(10.0,r/10.0) - 1.0);
B:= Ln(1/tanh(r/17.37) );
p:= sinh(B/(2*order));
p_2:= sqr(p);
for k:= 1 to order do
    begin
    ak[k]:= sin((2*k-1)*Pi/(2*order));
    bk[k]:= p_2 + sqr(sin(k*Pi/order));
    end;
g[1]:= 2*ak[1]/p;
for k:= 2 to order do
    begin
    g[k]:= 4*ak[k-1]*ak[k]/(bk[k-1]*g[k-1]);
    end;
if ( order mod 2 = 1 ) then
    g[order+1]:= 1.0 else g[order+1]:= sqr(1.0/tanh(0.25*B));
    gk_sum:= 0.0;
for k:= 1 to order do
    begin
    gk_sum:= gk_sum + g[k];
    end;
L_dB:= 4.343*gk_sum*sqrt(wp_lo*wp_hi)/( (wp_hi-wp_lo)*Qu );
T_do:= gk_sum/( 2*Pi*0.001*(wp_hi-wp_lo) );
phi:= 1/order*Arcsinh(1/eta);
if order mod 2 = 0 then Odd_order:= False else Odd_order:= True;
Const_2:= phi;
Const_1:= sinh(Const_2);
Const_2:= cosh(Const_2);

Peak:= 1.0;
if order mod 2 = 0 then
    begin
         for i:= 1 to order div 2 do
              begin
                   Pole_Re[i]:= -Const_1*sin(Pi*(2*i-1)/(2*order) );
                   Pole_Im[i]:=  Const_2*Cos(Pi*(2*i-1)/(2*order) );
                   Peak:= Peak * sqr( sqr(Pole_Re[i] ) + sqr(Pole_Im[i] ) );
              end;
         end
    else
         begin
              for i:= 1 to (order-1) div 2 do
                   begin
                   Pole_Re[i]:= -Const_1*sin(Pi*(2*i-1)/(2*order) );
                   Pole_Im[i]:=  Const_2*Cos(Pi*(2*i-1)/(2*order) );
                   Peak:= Peak * sqr( sqr(Pole_Re[i] ) + sqr(Pole_Im[i] ) );
              end;
         Pole_Re[i+1]:= -Const_1;
         Pole_Im[i+1]:= 0.0;
         Peak:= Peak * sqr( Const_1 );
    end;

w:= w_beg;
m:= 1;
wo_geom:= sqrt( wp_lo*wp_hi);
T_Lhs:= -1.0;
T_Rhs:= -1.0;
while w <= w_end do
    begin
    Time[m]:= 0.0;
    wp:= wo_geom/(wp_hi-wp_lo);
    wp:= wp*(w/wo_geom - wo_geom/w );
    Fr[m]:= w;
    if not Odd_order then
         begin
         Ang:= 0.0;
         for i:= 1 to order div 2 do
              begin
              t1:= wp - Pole_Im[i];
              t2:= wp + Pole_im[i];
              t3:= Pole_Re[i];
              Ang_1:= -Atan4(t1,t3);
              Time[m]:= Time[m] - sqr(Cos(Ang_1) ) /t3 ;
              Ang:= Ang + Ang_1;
              Ang_1:= -Atan4( t2,t3);
              Time[m]:= Time[m] - sqr(Cos(Ang_1) )/t3;
              Ang:= Ang + Ang_1;
              end;
         end;
    if Odd_order or (order = 1) then
    begin
         Ang:= 0.0;
         for i:= 1 to (order-1) div 2 do
              begin
              t1:= wp - Pole_Im[i];
              t2:= wp + Pole_im[i];
              t3:= Pole_Re[i];
              Ang_1:= -Atan4(t1,t3);
              Time[m]:= Time[m] - sqr(Cos(Ang_1) ) /t3 ;
              Ang:= Ang + Ang_1;
              Ang_1:= -Atan4( t2,t3);
              Time[m]:= Time[m] - sqr(Cos(Ang_1) )/t3;
              Ang:= Ang + Ang_1;
              end;
         t3:= Pole_Re[(order+1) div 2];
         Ang_1:= -Atan4(wp,t3);
         Time[m]:= Time[m] - sqr( Cos( Ang_1) ) /t3;
         Ang:= Ang + Ang_1;
         end;
Time[m]:= Time[m] * 1000 / ( Pi * (wp_hi-wp_lo) );
if w < wo_geom then
    begin
    if Time[m] > T_Lhs then T_Lhs:= Time[m];
    end;
if w > wo_geom then
    begin
    if Time[m] > T_Rhs then T_Rhs:= Time[m];
    end;
w:= w + w_inc;
m:= m + 1;
end;
w:= w_beg;
m:= 1;
clrscr;
while w < w_end do
    begin
    wp:= wo_geom/(wp_hi - wp_lo);
    wp:= abs( wp*(w/wo_geom - wo_geom/w ) );
    wp:= wp*1.0E10;
    wp:= Int(wp)*1.0E-10;
    if ( (wp > 1.0) and ( w < wp_lo) ) then
         begin
         mag[m]:= -10*Log10(1 + sqr(eta*Cosh(order*Arccosh(wp)))) - L_db*T_Lhs/T_do;
         end;
    if ( (wp > 1.0) and ( w > wp_hi) ) then
         begin
         mag[m]:= -10*Log10(1 + sqr(eta*Cosh(order*Arccosh(wp)))) - L_db*T_Rhs/T_do;
         end;
    if  wp <= 1.0 then
         begin
         mag[m]:= -10*Log10(1 + sqr(eta*Cos(order*Arccos(wp)))) - L_db*Time[m]/T_do;
         end;
if m <= 23 then
    begin
    gotoxy(1,m+1);
    write(Fr[m]:9:4,'  ',mag[m]:7:2);
    end;
if ( ( m > 23 ) and ( m <= 46) ) then
    begin
    gotoxy(20,m-22);
    write('  ',Fr[m]:9:4,'  ',mag[m]:7:2);
    end;
if ( ( m > 46 ) and ( m <= 69) ) then
    begin
    gotoxy(40,m-45);
    write('  ',Fr[m]:9:4,'  ',mag[m]:7:2);
    end;
if ( ( m > 69 ) and ( m <= 92) ) then
    begin
    gotoxy(60,m-68);
    write('  ',Fr[m]:9:4,' ',mag[m]:7:2);
    end;
w:= w + w_inc;
m:= m + 1;
end;
repeat until keypressed;
clrscr;
gotoxy(10,5);
write('Enter Title for Graph  ');
readln(Graph_Title);
for i:= 1 to m-1 do
    begin
    A[i,1]:= Fr[i];
    A[i,2]:= mag[i];
    end;
N:= m - 1;
assign(file_name,'PlotGr.Chn');
Chain(file_name);
readln(N);
clrscr;
gotoxy(10,5);
write('Enter Title for Graph for Group Delay ');
readln(Graph_Title);
for i:= 1 to m-1 do
    begin
    A[i,1]:= Fr[i];
    A[i,2]:= Time[i];
    end;
N:= m - 1;
assign(file_name,'PlotGr.Chn');
Chain(file_name);
Write_to_File( Fr,mag);
goto 500;
end.