#include
#include
#include
#define MAXN 10
void vrk4( double x0, double y0[], double h, int n, void (*fun)(double, double*, double*), double y1[] );
void Fpraw(double del,double y[], double f[]);
double Pierw(double M);
double siecz(double a,double b, double(*f)(double x), double eps);
double gam=1.4;
double M0=1.91;
double alfa0;
double costam;
double M;
int main()
{
alfa0=asin(1./M0);
double del, eps, a, b;
double ddel;
double y0[2], y1[2];
int n=2;

printf("Podaj poczatek i koniec przedzialu metody siecznych a i b:");
scanf("%lf%lf", &a, &b);
printf("\nPodaj krok metody siecznych eps: ");
scanf("%lf", &eps);
printf("\nPodaj krok całkowania h:");
scanf("%lf", &ddel);

// warunki poczatkowe
y0[0]= alfa0+sqrt((gam+1)/(gam-1))*atan(sqrt((gam-1)/
(gam+1))/tan(alfa0))-2.*atan(1.); //omega0
y0[1]=22.*atan(1.0)*4./180.; //teta0
del=40.*atan(1.0)*4./180.; //delta0=sigma
// calkowanie ukladu rownan
do
{
vrk4(del,y0,ddel,n,Fpraw,y1);
del=del+ddel;
printf("%lf %lf %lf\n",del,y1[0],y1[1]);
a=M-0.15;
b=M;
y0[0]=y1[0]; //omega
costam=y1[0];
y0[1]=y1[1]; //teta
M=siecz(a,b,Pierw,eps);
}
while(fabs(del-y1[1])>fabs(ddel));

printf("\nWynik: omega=%lf, teta=%lf\nKONIEC PROGRAMU\n", y1[0], 45*y1[1]/atan(1.));

system("PAUSE");
return 0;
}

void Fpraw(double del,double y[], double f[])
{
f[0]=(sin(y[1])*sin(del-y[1]))/
(tan(asin(1./M))*sin(del)*(1-(M*M*sin(del-y[1])*sin(del-y[1]))));
f[1]=-(sin(y[1])*cos(del-y[1]))/
(sin(del)*(1-(M*M*sin(del-y[1])*sin(del-y[1]))));
}


// wykonuje jeden krok calkowania
// wektorowego rownania rozniczkowego zwyczjanego:
//
// dY/dx = Fun(x,Y), Y(x0)=Y0
//
// metoda Rungego-Kutty IV-tego rzedu.
//
// Parametry formalne:
// x0 - wartosc startowa zm. niezaleznej
// y0 - wartosc startowa zm. zaleznej (tablica n-elementowa)
// h - krok calkowania
// n - liczba rownan
// fun(x,y,prawastr) - nazwa funkcji obliczajacej prawe strony
// y1 - obliczona wartosc zmiennej zaleznej w punkcie x0+h
// (tablica n-elementowa)

void vrk4( double x0, double y0[], double h, int n, void (*fun)(double, double*, double*), double y1[] )
{
int i;
double k1[MAXN], k2[MAXN], k3[MAXN], k4[MAXN];
double ytmp[MAXN];

fun( x0, y0, k1);
for ( i=0; i {
k1[i] *= h;
ytmp[i] = y0[i] + k1[i]/2.0;
}

fun( x0+h/2.0, ytmp, k2);
for ( i=0; i {
k2[i] *= h;
ytmp[i] = y0[i] + k2[i]/2.0;
}

fun( x0+h/2.0, ytmp, k3);
for ( i=0; i {
k3[i] *= h;
ytmp[i] = y0[i] + k3[i];
}

fun( x0+h, ytmp, k4);
for ( i=0; i k4[i] *= h;

for ( i=0; i y1[i] = y0[i] + (k1[i] + 2.*k2[i] + 2.*k3[i] + k4[i])/6.;

}
double Pierw(double M)
{
double pom;
pom=asin(1./M)+sqrt((gam+1)/(gam-1))*atan(sqrt((gam-1)/(gam+1))/tan(asin(1./M)))-atan(1.)*2.
-costam;
return pom;
}


double siecz(double a,double b, double(*f)(double x), double eps)
{
double x,x0,z;
x0=a;
do
{
x=a-f(a)*(b-a)/(f(b)-f(a));
if (f(a)*f(a)<0.0)
b=x;
else a=x;
z=fabs(x-x0);
x0=x;
}
while(z>eps);
return x;
}