Pole powierzchni xyz

Question

Witam, mam problem z napisaniem programu który obliczy pole powierzchni 400 punktów o współrzędnych xyz. Pole powierzchni ma być obliczane na zasadzie łączenia punktów w trójkąty, a następnie sumowaniu ich pól. Skorzystałem z programu Origin i okazało się, że wynik zwracany przez mój program jest błędny.

Wynik z Origin: 21712,737565358
Mój wynik: 22651.459434729  

Byłbym wdzięczny za wskazanie w czym może tkwić problem.

link do plik ze współrzędnymi - https://drive.google.com/file/d/1lJCaCc0hotaTqHJ7pp_khpwf2Jg9D3mj/view?usp=sharing

#include <iostream>
#include <fstream>
#include <iomanip>
#include <cmath>

using namespace std;

int main()  {

    ifstream wczytaj;
    wczytaj.open("calosc.dat");
    double A,B,C,A1,B1,C1,A2,B2,C2;
    double P;
    int n=400;
    double x[n],y[n],z[n];
    cout << setprecision(14);

    for(int i=0; i<n; i++)  {

        wczytaj >> x[i] >> y[i] >> z[i];
    }


    for(int i=0; i<2; i++)  { // Obliczenie 2 pierwszych trójkątów osobno

        if(i == 0){

            A=x[i+1]-x[i];
            B=y[i+1]-y[i];
            C=z[i+1]-z[i];

            A1=x[i+20]-x[i];
            B1=y[i+20]-y[i];
            C1=z[i+20]-z[i];

            A2=(B*C1)-(C*B1);
            B2=(A*C1)-(C*A1);
            C2=(A*B1)-(B*A1);

            P+=0.5*abs(sqrt(A2*A2+B2*B2+C2*C2));
        }

        if(i == 1){

            A=x[i+20]-x[i+19];
            B=y[i+20]-y[i+19];
            C=z[i+20]-z[i+19];

            A1=x[i]-x[i+20];
            B1=y[i]-y[i+20];
            C1=z[i]-z[i+20];

            A2=(B*C1)-(C*B1);
            B2=(A*C1)-(C*A1);
            C2=(A*B1)-(B*A1);

            P+=0.5*abs(sqrt(A2*A2+B2*B2+C2*C2));
        }
    }

    for(int i=2, k=0; i<399; i++,k++){

        if(i%2 == 0){

            A=x[i]-x[i-1];
            B=y[i]-y[i-1];
            C=z[i]-z[i-1];


            A1=x[i+19]-x[i-1];
            B1=y[i+19]-y[i-1];
            C1=z[i+19]-z[i-1];

            A2=(B*C1)-(C*B1);
            B2=(A*C1)-(C*A1);
            C2=(A*B1)-(B*A1);

            P+=0.5*abs(sqrt(A2*A2+B2*B2+C2*C2));
        }

        else if(i%2!=0){

            A=x[k+22]-x[k+21];
            B=y[k+22]-y[k+21];
            C=z[k+22]-z[k+21];

            A1=x[k]-x[k+21];
            B1=y[k]-y[k+21];
            C1=z[k]-z[k+21];

            A2=(B*C1)-(C*B1);
            B2=(A*C1)-(C*A1);
            C2=(A*B1)-(B*A1);

            P+=0.5*abs(sqrt(A2*A2+B2*B2+C2*C2));
        }
    }

    cout << "Pole powierzchni wynosi:  " << P << endl;

 return 0;
}