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

#define M2PI 6.28318530717959

int     i;
double  beta;
double  gama;
double  bmin;
double  bmax;
double  fdolni;
double  fhorni;
double *g;
double  R;
int     rada;
int     dolnip;


double dB2Np  (double x){return x*0.11512925464970229;}
double Np2dB  (double x){return x/0.11512925464970229;}
double MHz2pi (double x){return x*M2PI;}
double pi2MHz (double x){return x/M2PI;}
double acosh  (double x){return log(x+sqrt(x*x-1.));}

int main(int pocet,char *adresy[])
{
if (pocet<8 || pocet>9)
    {
    printf("\nProgram pro vypocet frekvencnich filtru - volne dle prilohy AR86 sepsal OK2PID\n");
    printf("\nParametry: H|D  C|B  R  dolni  horni  zvlneni  utlum  [rad] \n");
    printf("\nkde:");
    printf("\n   H = horni propust, D = dolni propust");
    printf("\n   C = Cebysevova aproximace, B = Butterwothova aproximace");
    printf("\n   R = vstupni impedance filtru (v ohmech)");
    printf("\n   dolni = hranice propustnosti filtru (v MHz)");
    printf("\n   horni = hranice nepropustnosti filtru (v MHz)");
    printf("\n   zvlneni = max. utlum (zvlneni) v propustne oblasti (v dB)");
    printf("\n   utlum = min. utlum v nepropustne oblasti (v dB)");
    printf("\n   rad = pozadovany rad filtru - pokud neni zadan, vypocita se\n\n");
    return EXIT_SUCCESS;
    }
switch (*adresy[1])
    {
    case 'd':
    case 'D':
	dolnip=1;
	break;
    case 'h':
    case 'H':
	dolnip=0;
	break;
    default:
	printf("%s? Neznamy typ propusti!",adresy[1]);
	return EXIT_SUCCESS;
    }
R     =      (atof(adresy[3]));
fdolni=MHz2pi(atof(adresy[4]));
fhorni=MHz2pi(atof(adresy[5]));
bmax  =dB2Np (atof(adresy[6]));
bmin  =dB2Np (atof(adresy[7]));
if (fdolni>=fhorni)
    {
    printf("\nDolni kmitocet nesmi byt vetsi nez horni!");
    return EXIT_SUCCESS;
    }
if (pocet==9)
    rada = atof(adresy[8]);
  else
    rada = 0;
 
switch (*adresy[2])    /* Typ filtru  */
    {
    case 'b':
    case 'B':
	printf("\nButterworthuv filtr");
	if (!rada)
	    rada = log ((exp(2.*bmin)-1.)/(exp(2.*bmax)-1.)) /
		   (2.* log(fhorni/fdolni));
	if (rada>8100)
	    {
	    printf("\nRad filtru vychazi prilis velky!\n");
	    return EXIT_SUCCESS;
	    }
	if (NULL==(g=malloc(sizeof(double)*(rada+3))))
	    {printf("Malo pameti!\n");return 1;}
	for (i=1,g[0]=0.;i<=rada;i++)
	    g[i]=2.*sin((2*i-1)*M_PI/(2.*rada));
	g[i]=1.;
	break;
    case 'c':
    case 'C':
	printf("\nCebysevuv filtr");
	if (!rada)
	    rada = acosh(sqrt((exp(2.*bmin-1.)/exp(2.*bmax)-1.))) /
		   acosh(fhorni/fdolni)+1;
	if (rada>8100)
	    {
	    printf("\n Rad filtru vychazi prilis velky!\n");
	    return EXIT_SUCCESS;
	    }
	if (NULL==(g=malloc(sizeof(double)*(rada+3))))
	    {printf("Malo pameti\n");return 1;}
	beta=log(1./tanh(bmax/2.));
	gama=sinh(beta/(2.*rada));
	g[0]=0.;
	g[1]=2.*sin(M_PI_2/rada) /gama;
	for (i=2;i<=rada;i++)
	    g[i]= 4. * sin((2*i-3)*M_PI/(2.*rada)) *
		       sin((2*i-1)*M_PI/(2.*rada)) /
		(g[i-1] * (gama*gama+sin((i-1)*M_PI/rada) *
				     sin((i-1)*M_PI/rada)));
	if (rada%2)
	    g[i]=1.;
	  else
	    g[i]=1./tanh(beta/4.) * 1./tanh(beta/4.);
	break;
    default:
	printf("%s? Neznamy typ filtru!",adresy[2]);
    }
printf (" ve tvaru PI");
if (dolnip)
    printf(" - dolni propust\n");
  else
    printf(" - horni propust\n");
printf("\nRad filtru:           %-2i"      ,rada);
printf("\nDolni kmitocet:       %5.3f MHz ",pi2MHz(fdolni));
printf("\nHorni kmitocet:       %5.3f MHz ",pi2MHz(fhorni));
printf("\nZvlneni:              %4.2f dB  ",Np2dB(bmax));
printf("\nUtlum:                %4.2f dB  ",Np2dB(bmin));
printf("\nVstupni impedance:   %3.0f ohmu, Vystupni impedance: %3.0f ohmu\n",
			R,R/g[rada+1]);

if (dolnip)
    for (i=1;i<=rada;i++)
	if (i%2)
	    printf("\nC%-2i = %6.0f pF",i,g[i]/R/fdolni*1e6);
	  else
	    printf("\nL%-2i = %6.3f uH",i,R*g[i]/fdolni);
  else
    for (i=1;i<=rada;i++)
	if (i%2)
	    printf("\nL%-2i = %6.3f uH",i,R/g[i]/fdolni);
	  else
	    printf("\nC%-2i = %6.0f pF",i,1e6/R/g[i]/fdolni);
printf("\n");
free (g);
return EXIT_SUCCESS;
}