Kolejne uruchomienia - różne wyniki

#include <iostream>
#include <fstream>
#include <iomanip>
#include <string>
#include"klasapanel.h"
using namespace std;
int main()
{
    //OTWARCIE PLIKÓW
    ifstream plik1;
    ofstream plik2;
    string odplik;
    string doplik;
    cout << "Podaj oznaczenie profilu\n";
    getline( cin, odplik );
    doplik = odplik + ".wir";
    odplik = odplik + ".stk";
    plik1.open( odplik.c_str() );
    //plik2.open(doplik.c_str());
    //DEKLARACJA I ZAPIS DANYCH GEOMETRYCZNYCH
    int ILE; //Ilość paneli modelujących profil
    plik1 >> ILE;
    float x1, x2, y1, y2; //Zmienne pomocnicze
    LinWir PWir[ ILE ]; //Tablica liniowych paneli wirowych
    //GEOMETRIA PANELI
    for( int i = 0; i < ILE; i++ )
    {
        plik1 >> x1 >> y1 >> x2 >> y2;
        PWir[ i ].PUNKTY( x1, y1, x2, y2 );
        cout << fixed << "\nPanel " << i << "\tx1= " << x1 << "\ty1= " << y1 << "\tx2= " << x2
        << "\ty2= " << y2 << "\tl= " << PWir[ i ].L();
    }
    plik1.close();
    //WARUNKI NAPŁYWU PŁYNU
    float alfa; //Geometryczny kąt natarcia
    float u, ux, uy; //Prędkość napływu płynu
    u = 10;
    cout << endl << "Podaj geometryczny kąt natarcia[deg]" << endl;
    cin >> alfa;
    alfa *= pi() / 360;
    ux = u * cos( alfa );
    uy = u * sin( alfa );
    //OBLICZENIA ROZKŁADU CYRKULACJI
    ///Macierz współczynników dla metody eliminacji Gaussa
    float ** A;
    A = new float *[ ILE + 1 ];
    for( int i = 0; i < ILE + 1; i++ )
    {
        A[ i ] = new float[ ILE + 1 ];
    }
    ///Zapisanie współczynników
    ///całki panelowe
    for( int i = 0; i < ILE; i++ )
    {
        for( int j = 0; j < ILE; j++ )
        {
            A[ i ][ j ] = PWir[ j ].CALPAN( PWir[ i ].XC(), PWir[ i ].YC() );
        }
    }
    ///przyczynki funkcji prądu i zera
    for( int i = 0; i < ILE; i++ )
    {
        A[ i ][ ILE ] = - 1;
        A[ ILE ][ i ] = 0;
    }
    ///warunek Kutty
    A[ ILE ][ 0 ] = 1, A[ ILE ][ ILE - 1 ] = 1;
    float X[ ILE + 1 ], Y[ ILE + 1 ];
    ///przyczynki przepływu prostego
    for( int i = 0; i < ILE; i++ )
    {
        Y[ i ] = - uy * PWir[ i ].XC() + ux * PWir[ i ].YC();
    }
    /*
    ////Zapis kontrolny*****************************
    for(int i=0; i<ILE+1; i++)
    {
    for(int j=0; j<ILE+1; j++)
    {
    plik2<<A[j][i]<<" ";
    }
    plik2<<Y[i];
    plik2<<endl;
    }
    plik2<<endl<<endl;
    for(int i=0; i<ILE; i++)
    {plik2<<PWir[i].XC()<<" ";}
    plik2<<endl;
    ////Koniec zapisu kontrolnego*****************************
    */
    ///Obliczenie rozkładu cyrkulacji metodą eliminacji Gaussa
    GAUSS( ILE + 1, A, X, Y );
    for( int i = 0; i < ILE; i++ )
    { plik2 << X[ i ] << " "; }
    x1 = 0;
    for( int i = 0; i < ILE; i++ )
    {
        PWir[ i ].gamma( X[ i ] );
        x1 = x1 + PWir[ i ].CYR();
    }
    //OBLICZENIE WSPÓŁCZYNNIKA SIŁY NOŚNEJ
    cout << "\nWspółczynnik siły nośnej to " << 2 * x1 /( u * PWir[ 0 ].XC() );
    ////////////////////////////////////////////////////////////////////////
    //SPRZĄTANIE PO PROGRAMIE
    for( int i = 0; i < ILE + 1; i++ ) delete[] A[ i ];
   
    { delete[] A; }
    return 0;
}

#include"MAT"
#include<iostream>
using namespace std;

//LINIA WIROWA - stały rozkład wirowości na linii prostej
class LinWir
{
public:
    //Konstruktor - goły
    LinWir();
    //Konstruktor - bez podania rozkładu cykrulacji
    LinWir( float x1, float y1, float x2, float y2 );
    //Konstruktor - z rozkładem cykrulacji
    LinWir( float x1, float y1, float x2, float y2, float cyr );
    //Zapis współrzędnych
    void PUNKTY( float x1, float y1, float x2, float y2 );
    //Zapis gęstości cyrkulacji
    void gamma( float cyr )
    {
        g = cyr;
        G = l * cyr;
    }
    //Całka do obliczenia funkcji prądu
    float CALPAN( float xr, float yr );
    //Wypis odciętej punktu kontrolnego
    float XC()
    { return xc; }
    //Wypis rzędnej punktu kontrolnego
    float YC()
    { return yc; }
    //Wypis gęstości cyrkulacji
    float cyr()
    { return g; }
    //Wypis długości profilu
    float L()
    { return l; }
    //Wypis cyrkulacji
    float CYR()
    { return G; }
private:
    //Współrzędne punktu początkowego, środkowego i końcowego
    float xp, yp, xc, yc, xk, yk;
    //Długość linii wirowej
    float l;
    //Liniowe natężenie cyrkulacji g=dG/dl
    float g;
    //Cyrkulacja
    float G;
};
///Definicja gołego konstruktora
LinWir::LinWir()
{ }
///Definicja konstruktora bez gęstości cyrkulacji
LinWir::LinWir( float x1, float y1, float x2, float y2 )
{
    xp = x1;
    xk = x2;
    yp = y1;
    yk = y2;
    xc = 0.5 *( xp + xk );
    yc = 0.5 *( yp + yk );
    l = pow( pow(( xk - xp ), 2 ) + pow(( yk - yp ), 2 ), 0.5 );
}
///Definicja konstruktora z gęstością cyrkulacji
LinWir::LinWir( float x1, float y1, float x2, float y2, float cyr )
{
    xp = x1;
    xk = x2;
    yp = y1;
    yk = y2;
    xc = 0.5 *( xp + xk );
    yc = 0.5 *( yp + yk );
    l = pow( pow(( xk - xp ), 2 ) + pow(( yk - yp ), 2 ), 0.5 );
    g = cyr;
    G = g * l;
}
///Definicja funkcji PUNKTY
void LinWir::PUNKTY( float x1, float y1, float x2, float y2 )
{
    xp = x1;
    xk = x2;
    yp = y1;
    yk = y2;
    xc = 0.5 *( xp + xk );
    yc = 0.5 *( yp + yk );
    l = pow( pow(( xk - xp ), 2 ) + pow(( yk - yp ), 2 ), 0.5 );
}
///Definicja funkcji CALPAN
float LinWir::CALPAN( float xr, float yr )
{
    int dok = 11;
    float xa[ dok + 1 ], ya[ dok + 1 ], F[ dok + 1 ], dx2, dy2;
    float sinfi =( yk - yp ) / l;
    float cosfi =( xk - xp ) / l;
    float lx[ dok + 1 ];
    for( int i = 0; i < dok + 1; i++ )
    {
        lx[ i ] = i * l /( dok );
        xa[ i ] = xp + cosfi * lx[ i ];
        ya[ i ] = yp + sinfi * lx[ i ];
        dx2 = pow( xr - xa[ i ], 2 );
        dy2 = pow( yr - ya[ i ], 2 );
        F[ i ] = log( pow( dx2 + dy2, 0.5 ) );
    }
    return - CALKA( dok + 1, lx, F ) /( 2 * pi() );
}