Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- #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 300 //в Td
- #define d 0.3918e-4
- #define T0 0.26 //500K 1K = 1.3806*10^-16 erg
- #define eps0 1.6e-12 //перевод энергии
- #define B0 2.4668e-4
- #define Tk 0.3
- #define hw 0.29299
- #define N 100001 //кол-во шагов по времени
- #define M 200 //кол-во шагов по энергии
- #define t 1e-14 //шаг по времени
- #define h 0.1 //шаг по энергии эВ
- double T[100][2]; // транспортное сечение
- double R[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");
- for (int i = 0; i<100; i++) {
- for (int j = 0; j<2; j++) {
- T[i][j] = 0;
- R[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]);
- }
- }
- fclose(transp);
- fclose(rot);
- }
- 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;
- }
- }
- return sigma;
- }
- void main() {
- FILE *FREE; FREE = fopen("C:\\FREE.txt", "w");
- double A, B, C, D, E; // константы в изначальной разностной схеме
- 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;
- double Fp[M], Fn[M]; // Fp - предыдуший шаг по времени, Fn - следующий шаг по времени
- /* double **F = new double*[3]; // 3 столбцов (время), M строк (энергия)
- for (int i = 0; i < 3; i++) { // массив для шагов по времени 10,50,100,500,1000.........
- F[i] = new double[M];
- }*/
- double F[6][M];
- for (int j = 0; j < M; j++) { //начальное условие
- Fp[j] = exp(-h*j / T0);
- }
- // Fp[M - 1]; // ГУ справа
- // Fp[0] = Fp[1]; // ГУ слева
- for (int i = 0; i < N; 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'));
- Fn[j] = (1 / (A*sqrt(h*j))) * (Fp[j + 1] * t* (C1 / (h*h) + C3 / h) +
- +Fp[j] * (A*sqrt(h*j) - 2 * t*C1 / (h*h) - t*C3 / h + t*C4) + Fp[j - 1] * (t*C1 / (h*h)));
- }
- Fn[M - 1] = 0;
- Fn[0] = Fn[1];
- for (int j = 0; j < M; j++) {
- Fp[j] = Fn[j];
- if (i == 10) F[0][j] = Fn[j];
- if (i == 50) F[1][j] = Fn[j];
- if (i == 100) F[2][j] = Fn[j];
- if (i == 1000) F[3][j] = Fn[j];
- if (i == 10000) F[4][j] = Fn[j];
- if (i == 100000) F[5][j] = Fn[j];
- // printf("%e\n", F[0][j]);
- // Sleep(15);
- }
- // printf("\n");
- }
- for (int j = 0; j < M; j++) {
- for (int i = 0; i < 6; i = i++) {
- fprintf(FREE, "%e\t", F[i][j]);
- }
- fprintf(FREE, "\n");
- }
- fclose(FREE);
- }
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement
