VALEK

#define _CRT_SECURE_NO_WARNINGS
#include <iostream>
#include <stdio.h>
#include <cmath>
#include <windows.h>
#define _USE_MATH_DEFINES
using namespace std;

#define m 9.1094e-28
#define MN2 46.4968e-24
#define n0 1e16
#define q 4.8032e-10
#define En 70           //в Td
#define d 0.3918e-4
#define T0 0.26 //300K 1K = 1.3806*10^-16 erg
#define eps0 1.6e-12 //перевод энергии
#define B0 2.4668e-4
#define hw 0.29299 // для вычисление колебательного уровня, гармоническая контанства в Эв
#define Tk 0.3  //3000K


#define N 10000 //кол-во шагов по времени
#define M 150 //кол-во шагов по энергии
#define t 1e-12 //шаг по времени
#define h 0.1 //шаг по энергии  эВ



double T[100][2];                // объявил глобальный массив
double R[100][2];
double VL[100][2];
double e;                        // передаваемое значение энергии
char type;                      // тип сечения
int a;

void init() {

    FILE *transp; transp = fopen("C:\\Users\\superior\\Documents\\cross_section\\transfer.txt", "r");
    FILE *rot; rot = fopen("C:\\Users\\superior\\Documents\\cross_section\\rotat.txt", "r");
    FILE *vl; vl = fopen("C:\\Users\\superior\\Documents\\cross_section\\vl.txt", "r");

    for (int i = 0; i<100; i++) {
        for (int j = 0; j<2; j++) {
            T[i][j] = 0;
            R[i][j] = 0;
            VL[i][j] = 0;
        }
    }

    for (int i = 0; i<100; i++) {
        for (int j = 0; j<2; j++) {
            fscanf(transp, "%lf", &T[i][j]);
            fscanf(rot, "%lf", &R[i][j]);
            fscanf(vl, "%lf", &VL[i][j]);
        }
    }

    fclose(transp);
    fclose(rot);
    fclose(vl);
}

double cs(double e, char type) {
    double sigma;

    if (type == 'r') {
        for (int i = 0; i < 45; i++) {
            if (R[i][0] >= e) {
                sigma = (R[i - 1][1] + (R[i][1] - R[i - 1][1])*(e - R[i - 1][0]) / (R[i][0] - R[i - 1][0]))*1e-16;
                break;
            }
        }

        if (R[44][0] < e) {
            sigma = (R[43][1] + (R[44][1] - R[43][1])*(e - R[43][0]) / (R[44][0] - R[43][0]))*1e-16;
        }
        if (R[0][0] > e) {
            sigma = (R[0][1] + (R[1][1] - R[0][1])*(e - R[0][0]) / (R[1][0] - R[0][0]))*1e-16;
        }
    }

    else if (type == 't') {
        for (int i = 0; i < 64; i++) {
            if (T[i][0] >= e) {
                sigma = (T[i - 1][1] + (T[i][1] - T[i - 1][1])*(e - T[i - 1][0]) / (T[i][0] - T[i - 1][0]))*1e-16;
                break;
            }

        }
        if (T[63][0] < e) {
            sigma = (T[62][1] + (T[63][1] - T[62][1])*(e - T[62][0]) / (T[63][0] - T[62][0]))*1e-16;
        }
        if (T[0][0] > e) {
            sigma = (T[0][1] + (T[1][1] - T[0][1])*(e - T[0][0]) / (T[1][0] - T[0][0]))*1e-16;
        }
    }
    else if (type == 'l') {
            for (int i = 0; i < 68; i++) {
                if (VL[i][0] >= e) {
                    sigma = (VL[i - 1][1] + (VL[i][1] - VL[i - 1][1])*(e - VL[i - 1][0]) / (VL[i][0] - VL[i - 1][0]))*1e-16;
                    break;
                }

            }
            if (VL[67][0] < e) {
                sigma = (VL[66][1] + (VL[67][1] - VL[66][1])*(e - VL[66][0]) / (VL[67][0] - VL[66][0]))*1e-16;
            }
            if (VL[0][0] > e) {
                sigma = (VL[0][1] + (VL[1][1] - VL[0][1])*(e - VL[0][0]) / (VL[1][0] - VL[0][0]))*1e-16;
            }

        }

        return sigma;
    }




void main() {
    FILE *FREE; FREE = fopen("C:\\FREE.txt", "w");

    double A,B,C,D,E, G;            // константы в изначальной разностной схеме
    double H;                   // член с электронным возбуждением
    double C1,C2,C3,C4;     // контанты в конечной разностной схеме
    double I=0;             // нормировочный интеграл

    init();

    A = (1 / n0)*sqrt(m / 2)*sqrt(eps0);
    B = ((q*En*1e-19)/3)*((q*En*1e-19)/3)/(3*eps0) ;
    C = d*T0*eps0;
    D = d*eps0;
    E = 4 * B0*eps0;
    G = exp(-hw / (2 * Tk))*eps0;


    double **F = new double*[N]; //N строк (время), M столбцов (энергия)
    for (int i = 0; i < N; i++) {
        F[i] = new double[M];
    }
    double **F0 = new double*[N]; //N строк (время), M столбцов (энергия)
    for (int i = 0; i < N; i++) {   // нормированная функция распределения
        F0[i] = new double[M];
    }


    for (int j = 0; j < M; j++) {   //начальное условие
        F[0][j] = exp(-h*j / T0);
    }
    for (int i = 0; i < N; i++) {   // граничное условие справа
        F[i][M-1] = 0;
    }

    int Evl;

    Evl = round(7.65/h);

    for (int i = 0; i < N-1 ; i++) {
        for (int j = 1; j < M-1; j++) {


            C1 = B*(h*j) / cs(h*j, 't') + C*(h*j)*(h*j)*cs(h*j, 't');
            C2 = D*(h*j)*(h*j)*cs(h*j, 't') + E*(h*j)*cs(h*j, 'r');   // !!!!!!!!!!!!!! тут второе сечение должно быть вращательным
            C3 = B*(cs(h*j, 't') - j*(cs(h*j + 0.5*h, 't') - cs(h*j - 0.5*h, 't'))) / (cs(h*j, 't')*cs(h*j, 't')) +
                + 2 * C*(h*j)*cs(h*j, 't') + C*(h*j)*(j)*(cs(h*j + 0.5*h, 't') - cs(h*j - 0.5*h, 't')) + C2;
            C4 = 2 * D*(h*j)*cs(h*j, 't') + D*(h*j)*(j)*(cs(h*j + 0.5*h, 't') - cs(h*j - 0.5*h, 't')) + E*cs(h*j, 'r') +
                + E*(h*j)*(cs(h*j + 0.5*h, 'r') - cs(h*j - 0.5*h, 'r')); // внимание - два последих сечения - вращательные!!!

            if (j + Evl < M) {
                H = (1 / A)*(1 / sqrt(h*j))*t*G*(h*j + Evl*h)*cs(h*j + h*Evl, 'l');
            }
            else H = 0;

            F[i + 1][j] = (1 / (A*sqrt(h*j))) * (F[i][j + 1] * t* (C1 / (h*h) + C3 / h) +
                +F[i][j] * (A*sqrt(h*j) - 2 * t*C1 / (h*h) - t*C3 / h + t*C4 - t*G*(h*j)*cs(h*j, 'l')) + F[i][j - 1] * (t*C1 / (h*h))) + H*F[i][j + Evl];

//          printf("%e\n", H);

//          Sleep(100);



        }
//      printf("\n");
        F[i+1][0] = F[i+1][1];                                  //ГУ слева
        F[0][0] = F[0][1];
    }
// нормировка
    for (int i = 0; i < N; i++) {
        for (int j = 1; j < M; j++) {
            I = I + (F[i][j - 1] + F[i][j])*h / 2;
        }
        for (int L = 0; L < M; L++) {
            F0[i][L] = F[i][L] / I;
        }
//      printf("%e\n", Evl);
//      Sleep(1500);

        I = 0;
    }

//
    for (int j = 0; j < M; j++) {

        for (int i = 1; i < N; i++) {

//      for (int i = 0; i < 5; i++) {
                fprintf(FREE, "%e\t", F0[i][j]);
            }

        }
        fprintf(FREE, "\n");
    }



/*
типы сечений:
Транспортное = t;
Вращательное = r;
электронное = l;
колебательное = w;
*/