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

#include "../src/nrutil.c"
#include "../src/nrutil.h"
#include "../src/gaussj.c"

int main(void)
{
    int i, j, k, n=200;
    float q, q1=0.2, q2=5.0, o;
    float **a, **b, **c, **tmpa, **tmpb;
    a=matrix(1,5,1,5);
    tmpa=matrix(1,5,1,5);
    b=matrix(1,5,1,1);
    tmpb=matrix(1,5,1,1);
    c=matrix(1,5,1,1);


    FILE *file;
    file=fopen("dane.dat","w");

    if(!file)
        printf("\nNie mozna otworzyc pliku!");

    for(i=1;i<=n;++i)
    {
        q=q1+((q2-q1)/n*(i-1));

        a[1][1]=q*0.002;
        a[1][2]=1;
        a[1][3]=6;
        a[1][4]=9;
        a[1][5]=10;
        a[2][1]=0.002;
        a[2][2]=1;
        a[2][3]=6;
        a[2][4]=9;
        a[2][5]=10;
        a[3][1]=1;
        a[3][2]=6;
        a[3][3]=6;
        a[3][4]=8;
        a[3][5]=6;
        a[4][1]=5;
        a[4][2]=9;
        a[4][3]=10;
        a[4][4]=7;
        a[4][5]=10;
        a[5][1]=3;
        a[5][2]=4;
        a[5][3]=9;
        a[5][4]=7;
        a[5][5]=9;

        for(j=1;j<=5;++j)
            for(k=1;k<=5;++k)
                tmpa[j][k]=a[j][k];

        b[1][1]=10;
        b[2][1]=2;
        b[3][1]=9;
        b[4][1]=9;
        b[5][1]=3;

        for(j=1;j<=5;++j)
            for(k=1;k<=1;++k)
                tmpb[j][k]=b[j][k];

        gaussj(a,5,b,1);

        for(j=1;j<=5;++j)
        {
            c[j][1]=0.0;
            for(k=1;k<=5;++k)
                c[j][1]+=tmpa[j][k]*b[k][1];
        }

        float tmpo=0;
        for(j=1;j<=5;++j)
            tmpo+=pow((c[j][1]-tmpb[j][1]),2);

        o=0.2*sqrt(tmpo);
        fprintf(file,"q=%f o=%f\n",q, o);
    }

    free_matrix(a,1,5,1,5);
    free_matrix(b,1,5,1,1);
    free_matrix(c,1,5,1,1);
    free_matrix(tmpa,1,5,1,5);
    free_matrix(tmpb,1,5,1,1);
    fclose(file);

    return 0;
}