1.cpp

#include <iostream>
//#include "math.h"
#include <vector>
#include <fstream>
#include <ctime>

using namespace std;

void readFile();    // Read data from file
void calculate();
int sign(double);   // signed function

vector<double> m, Eg, Ec, dSO, En;

int main() {
// read data from file

    readFile();

// reference point of time

    clock_t begin = clock();

// the function which shoud be improved

    calculate();

// end of time
clock_t end = clock();

// spent time (sec)
double ec_time = double(end - begin) / CLOCKS_PER_SEC;

cout << "time: " << ec_time << "\n";

cout << En[0] << "\n"
     << En[1] << "\n"
     << En[2] << "\n"
     << En[3] << "\n"
     << En[4] << "\n";
}

int sign (double Value){
    int sValue=-1;
    if (Value > 0) sValue=1;
    return sValue;
}

void readFile() {

    ifstream infile;
            infile.open ("1.txt");
    double a, b, c, d;
    m.clear (); Eg.clear (); Ec.clear (); dSO.clear ();
    while (infile >> a >> b >> c >> d)
    {
     m.push_back(a);
     Ec.push_back(b);
     Eg.push_back(c);
     dSO.push_back(d/3+c);
    }
    infile.close();
}

void calculate() {
    double pi=3.141592653589793;
    double nm=1e-9;
    double h=6.626068e-34/(2*pi);
    double e=1.60217646e-19;
    double dz=0.0212306; // the step between points
    double D_const = dz*dz*nm*nm/(h*h); // some const
    int setWFs=50; // the number of required zero values
    double x1, x2; // the values of "new_val" between different E
    double dE=e*1e-5; // the step of E

const size_t num_poi=m.size (); // the number of z-points
vector<double> m_E(num_poi); // the vector of m which dependent from E
double E=3.74983e-19; // min(Ec); the minimum value of Ec

//calculation of vector m(E) at start point min(Ec)
for (size_t i = 0; i < num_poi; ++i)  m_E[i]=m[i]*(1-(Ec[i]-E)/dSO[i]);

// initial values of new_val at start z-points
double new_val=0, prev_val = 1.0, prev_2_val = 0.0;

// tha calculation of new_val at last z-point
for (size_t i = 2; i < num_poi; ++i) {
    new_val =
         (2*D_const*(Ec[i]-E)+
          2/(m_E[i]+m_E[i-1])+
          2/(m_E[i-1]+m_E[i-2]))*(m_E[i]+m_E[i-1])*prev_val/2+
            (m_E[i]+m_E[i-1])*(-prev_2_val)/((m_E[i-1]+m_E[i-2]));
                    prev_2_val = prev_val;
                    prev_val = new_val;
}

// next step of E
E+=dE;

// clear information about found zero-points
En.clear ();

// similar calculation until En.size!=setWFs
do {
x1=new_val;

    for (size_t i = 0; i < num_poi; ++i)  m_E[i]=m[i]*(1-(Ec[i]-E)/dSO[i]);

  prev_val = 1.0, prev_2_val = 0.0;
                  for (size_t i = 2; i !=num_poi; ++i) {
  new_val =
            (2*D_const*(Ec[i]-E)+
             2/(m_E[i]+m_E[i-1])+
             2/(m_E[i-1]+m_E[i-2]))*(m_E[i]+m_E[i-1])*prev_val/2+
             (m_E[i]+m_E[i-1])*(-prev_2_val)/((m_E[i-1]+m_E[i-2]));
                          prev_2_val = prev_val;
                          prev_val = new_val;
             }

x2=new_val;
// if sign(x) isn't changed -> next E = E + dE;
        if (sign(x1)==sign(x2)) E+=dE;
// if sign(x) is chenged -> write value of E and -> next E = E +dE
        else {

            En.push_back (E-dE-x1*dE/(x2-x1));
            E+=dE;
            cout << "Found\t" << En.size () << "\tof\t" << setWFs << "\n";

        }

} while (En.size()<setWFs);

}