API tools faq
paste
phystota

FREEv3.0

phystota
Apr 13th, 2019
38
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
C 26.39 KB | None | 0 0
raw download clone embed print report
  1. #define _CRT_SECURE_NO_WARNINGS
  2. #include <iostream>
  3. #include <stdio.h>
  4. #include <cmath>
  5. #include <windows.h>
  6. #define _USE_MATH_DEFINES
  7. using namespace std;
  8.  
  9. #define m 9.1094e-28
  10. #define MN2 46.4968e-24
  11. #define n0 1e16
  12. #define ne 1e12
  13. #define q 4.8032e-10
  14. #define En 0          //в Td    
  15. #define d 3.91665e-5
  16. #define T0 0.02586 //300K 1K = 1.3806*10^-16 erg
  17. #define Te 1.5
  18. #define eps0 1.6e-12 //перевод энергии
  19. #define B0 2.4668e-4
  20. #define Tk 0.25862 //3000K
  21. #define hw 0.29299
  22. #define hw0 0.00178
  23.  
  24. #define logL 10
  25. #define PI 3.14159265
  26.  
  27.  
  28. #define N 100001 //кол-во шагов по времени
  29. #define M 400 //кол-во шагов по энергии
  30. #define t 1e-16//шаг по времени
  31. #define h 0.1 //шаг по энергии  эВ
  32. #define z 8 //количество нужных шагов по энергии
  33.  
  34. double T[100][2];               // транспортное сечение
  35. double R[100][2];               // вращательное сечение
  36. double VL[100][2];              // электронное возбуждение
  37. double VW[100][2];              // полная колебаловка
  38.  
  39. double W1[100][2], W2[100][2], W3[100][2], W4[100][2], W5[100][2], W6[100][2], W7[100][2], W8[100][2];      // колебаловка поуровневая
  40.  
  41. double CS_Ex_1[100][2], CS_Ex_2[100][2], CS_Ex_3[100][2], CS_Ex_4[100][2], CS_Ex_5[100][2],
  42. CS_Ex_6[100][2], CS_Ex_7[100][2], CS_Ex_8[100][2], CS_Ex_9[100][2], CS_Ex_10[100][2], CS_Ex_11[100][2]; // электронное возбуждеиние
  43.  
  44. double I11, I21, I13, I14, I15, I16, I17, I18 = 0; // интегралы при решении уравнения для ЗС колебательной энергии
  45.  
  46. double e;                        // передаваемое значение энергии  
  47. char type;                      // тип сечения
  48. int a;
  49.  
  50. void init() {
  51.  
  52.     FILE *transp; transp = fopen("C:\\Users\\superior\\Documents\\cross_section\\transfer.txt", "r");
  53.     FILE *rot; rot = fopen("C:\\Users\\superior\\Documents\\cross_section\\rotat.txt", "r");
  54.     FILE *vl; vl = fopen("C:\\Users\\superior\\Documents\\cross_section\\vl.txt", "r");
  55.     FILE *vw; vw = fopen("C:\\Users\\superior\\Documents\\cross_section\\vw.txt", "r");
  56.  
  57.     FILE *w1; w1 = fopen("C:\\Users\\superior\\Documents\\cross_section\\w1.txt", "r");
  58.     FILE *w2; w2 = fopen("C:\\Users\\superior\\Documents\\cross_section\\w2.txt", "r");
  59.     FILE *w3; w3 = fopen("C:\\Users\\superior\\Documents\\cross_section\\w3.txt", "r");
  60.     FILE *w4; w4 = fopen("C:\\Users\\superior\\Documents\\cross_section\\w4.txt", "r");
  61.     FILE *w5; w5 = fopen("C:\\Users\\superior\\Documents\\cross_section\\w5.txt", "r");
  62.     FILE *w6; w6 = fopen("C:\\Users\\superior\\Documents\\cross_section\\w6.txt", "r");
  63.     FILE *w7; w7 = fopen("C:\\Users\\superior\\Documents\\cross_section\\w7.txt", "r");
  64.     FILE *w8; w8 = fopen("C:\\Users\\superior\\Documents\\cross_section\\w8.txt", "r");
  65.  
  66.     FILE *CS_Ex_1_file; CS_Ex_1_file = fopen("C:\\cross_section\\ex\\E1.txt", "r");
  67.     FILE *CS_Ex_2_file; CS_Ex_2_file = fopen("C:\\cross_section\\ex\\E2.txt", "r");
  68.     FILE *CS_Ex_3_file; CS_Ex_3_file = fopen("C:\\cross_section\\ex\\E3.txt", "r");
  69.     FILE *CS_Ex_4_file; CS_Ex_4_file = fopen("C:\\cross_section\\ex\\E4.txt", "r");
  70.     FILE *CS_Ex_5_file; CS_Ex_5_file = fopen("C:\\cross_section\\ex\\E5.txt", "r");
  71.     FILE *CS_Ex_6_file; CS_Ex_6_file = fopen("C:\\cross_section\\ex\\E6.txt", "r");
  72.     FILE *CS_Ex_7_file; CS_Ex_7_file = fopen("C:\\cross_section\\ex\\E7.txt", "r");
  73.     FILE *CS_Ex_8_file; CS_Ex_8_file = fopen("C:\\cross_section\\ex\\E8.txt", "r");
  74.     FILE *CS_Ex_9_file; CS_Ex_9_file = fopen("C:\\cross_section\\ex\\E9.txt", "r");
  75.     FILE *CS_Ex_10_file; CS_Ex_10_file = fopen("C:\\cross_section\\ex\\E10.txt", "r");
  76.     FILE *CS_Ex_11_file; CS_Ex_11_file = fopen("C:\\cross_section\\ex\\E11.txt", "r");
  77.  
  78.  
  79.     for (int i = 0; i < 100; i++) {
  80.         for (int j = 0; j < 2; j++) {
  81.             T[i][j] = 0; R[i][j] = 0; VL[i][j] = 0; VW[i][j] = 0;
  82.  
  83.             W1[i][j] = 0;W2[i][j] = 0; W3[i][j] = 0; W4[i][j] = 0; W5[i][j] = 0; W6[i][j] = 0; W7[i][j] = 0; W8[i][j] = 0;
  84.  
  85.             CS_Ex_1[i][j] = 0; CS_Ex_2[i][j] = 0; CS_Ex_3[i][j] = 0; CS_Ex_4[i][j] = 0; CS_Ex_5[i][j] = 0; CS_Ex_6[i][j] = 0; CS_Ex_7[i][j] = 0; CS_Ex_8[i][j] = 0; CS_Ex_9[i][j] = 0; CS_Ex_10[i][j] = 0; CS_Ex_11[i][j] = 0;
  86.         }
  87.     }
  88.  
  89.  
  90.     for (int i = 0; i<100; i++) {
  91.         for (int j = 0; j<2; j++) {
  92.             fscanf(transp, "%lf", &T[i][j]);
  93.             fscanf(rot, "%lf", &R[i][j]);
  94.             fscanf(vl, "%lf", &VL[i][j]);
  95.             fscanf(vw, "%lf", &VW[i][j]);
  96.            
  97.             fscanf(w1, "%lf", &W1[i][j]);
  98.             fscanf(w2, "%lf", &W2[i][j]);
  99.             fscanf(w3, "%lf", &W3[i][j]);
  100.             fscanf(w4, "%lf", &W4[i][j]);
  101.             fscanf(w5, "%lf", &W5[i][j]);
  102.             fscanf(w6, "%lf", &W6[i][j]);
  103.             fscanf(w7, "%lf", &W7[i][j]);
  104.             fscanf(w8, "%lf", &W8[i][j]);
  105.  
  106.             fscanf(CS_Ex_1_file, "%lf", &CS_Ex_1[i][j]);
  107.             fscanf(CS_Ex_2_file, "%lf", &CS_Ex_2[i][j]);
  108.             fscanf(CS_Ex_3_file, "%lf", &CS_Ex_3[i][j]);
  109.             fscanf(CS_Ex_4_file, "%lf", &CS_Ex_4[i][j]);
  110.             fscanf(CS_Ex_5_file, "%lf", &CS_Ex_5[i][j]);
  111.             fscanf(CS_Ex_6_file, "%lf", &CS_Ex_6[i][j]);
  112.             fscanf(CS_Ex_7_file, "%lf", &CS_Ex_7[i][j]);
  113.             fscanf(CS_Ex_8_file, "%lf", &CS_Ex_8[i][j]);
  114.             fscanf(CS_Ex_9_file, "%lf", &CS_Ex_9[i][j]);
  115.             fscanf(CS_Ex_10_file, "%lf", &CS_Ex_10[i][j]);
  116.             fscanf(CS_Ex_11_file, "%lf", &CS_Ex_11[i][j]);
  117.         }
  118.     }
  119.  
  120.     fclose(transp); fclose(rot);  fclose(vl); fclose(vw);
  121.  
  122.     fclose(w1); fclose(w2); fclose(w3); fclose(w4); fclose(w5); fclose(w6); fclose(w7); fclose(w8);
  123.  
  124.     fclose(CS_Ex_1_file); fclose(CS_Ex_2_file); fclose(CS_Ex_3_file); fclose(CS_Ex_4_file); fclose(CS_Ex_5_file); fclose(CS_Ex_6_file); fclose(CS_Ex_7_file); fclose(CS_Ex_8_file); fclose(CS_Ex_9_file); fclose(CS_Ex_10_file); fclose(CS_Ex_11_file);
  125.  
  126. }
  127.  
  128.  
  129. double cs(double e, string type) {
  130.     double sigma;
  131.  
  132.     if (type == "r") {
  133.         for (int i = 0; i < 26; i++) {
  134.             if (R[i][0] == e) {
  135.                 sigma = R[i][1];
  136.                 return sigma;
  137.             }
  138.  
  139.             if (R[i][0] > e) {
  140.                 sigma = (R[i - 1][1] + (R[i][1] - R[i - 1][1])*(e - R[i - 1][0]) / (R[i][0] - R[i - 1][0]));
  141.                 return sigma;
  142.             }
  143.  
  144.         }
  145.         if (R[25][0] < e) {
  146.             sigma = 0;
  147.         }
  148.         if (R[0][0] > e) {
  149.             sigma = 0;
  150.         }
  151.     }
  152.  
  153.     else if (type == "t") {
  154.         for (int i = 0; i < 54; i++) {
  155.             if (T[i][0] == e) {
  156.                 sigma = T[i][1];
  157.                 return sigma;
  158.             }
  159.  
  160.             if (T[i][0] > e) {
  161.                 sigma = (T[i - 1][1] + (T[i][1] - T[i - 1][1])*(e - T[i - 1][0]) / (T[i][0] - T[i - 1][0]));
  162.                 return sigma;
  163.             }
  164.  
  165.         }
  166.         if (T[53][0] < e) {
  167.             sigma = 0;
  168.         }
  169.         if (T[0][0] > e) {
  170.             sigma = 0;
  171.         }
  172.     }
  173.     else if (type == "l") {
  174.         for (int i = 0; i < 68; i++) {
  175.             if (VL[i][0] > e) {
  176.                 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;
  177.                 break;
  178.             }
  179.             if (VL[i][0] == e) {
  180.                 sigma = VL[i][1] * 1e-16;
  181.                 break;
  182.             }
  183.  
  184.         }
  185.         if (VL[66][0] < e) {
  186.             sigma = (VL[65][1] + (VL[66][1] - VL[65][1])*(e - VL[65][0]) / (VL[66][0] - VL[65][0]))*1e-16;
  187.         }
  188.  
  189.     }
  190.     else if (type == "w") {
  191.         for (int i = 0; i < 28; i++) {
  192.             if (VW[i][0] == e) {
  193.                 sigma = VW[i][1];
  194.                 return sigma;
  195.             }
  196.  
  197.             if (R[i][0] > e) {
  198.                 sigma = (R[i - 1][1] + (R[i][1] - R[i - 1][1])*(e - R[i - 1][0]) / (R[i][0] - R[i - 1][0]));
  199.                 return sigma;
  200.             }
  201.  
  202.         }
  203.         if (R[25][0] < e) {
  204.             sigma = 0;
  205.         }
  206.         if (R[0][0] > e) {
  207.             sigma = 0;
  208.         }
  209.     }
  210.     else if (type == "w1") {
  211.         for (int i = 0; i < 49; i++) {
  212.             if (W1[i][0] == e) {
  213.                 sigma = W1[i][1];
  214.                 return sigma;
  215.             }
  216.  
  217.             if (W1[i][0] > e) {
  218.                 sigma = (W1[i - 1][1] + (W1[i][1] - W1[i - 1][1])*(e - W1[i - 1][0]) / (W1[i][0] - W1[i - 1][0]));
  219.                 return sigma;
  220.             }
  221.  
  222.         }
  223.         if (W1[48][0] < e) {
  224.             sigma = 0;
  225.             return sigma;
  226.         }
  227.         if (W1[0][0] > e) {
  228.             sigma = 0;
  229.             return sigma;
  230.         }
  231.     }
  232.     else if (type == "w2") {
  233.         for (int i = 0; i < 24; i++) {
  234.             if (W2[i][0] == e) {
  235.                 sigma = W2[i][1];
  236.                 return sigma;
  237.             }
  238.  
  239.             if (W2[i][0] > e) {
  240.                 sigma = (W2[i - 1][1] + (W2[i][1] - W2[i - 1][1])*(e - W2[i - 1][0]) / (W2[i][0] - W2[i - 1][0]));
  241.                 return sigma;
  242.             }
  243.  
  244.         }
  245.         if (W2[23][0] < e) {
  246.             sigma = 0;
  247.             return sigma;
  248.         }
  249.         if (W2[0][0] > e) {
  250.             sigma = 0;
  251.             return sigma;
  252.         }
  253.     }
  254.     else if (type == "w3") {
  255.         for (int i = 0; i < 20; i++) {
  256.             if (W3[i][0] == e) {
  257.                 sigma = W3[i][1];
  258.                 return sigma;
  259.             }
  260.  
  261.             if (W3[i][0] > e) {
  262.                 sigma = (W3[i - 1][1] + (W3[i][1] - W3[i - 1][1])*(e - W3[i - 1][0]) / (W3[i][0] - W3[i - 1][0]));
  263.                 return sigma;
  264.             }
  265.  
  266.         }
  267.         if (W3[19][0] < e) {
  268.             sigma = 0;
  269.             return sigma;
  270.         }
  271.         if (W3[0][0] > e) {
  272.             sigma = 0;
  273.             return sigma;
  274.         }
  275.     }
  276.     else if (type == "w4") {
  277.         for (int i = 0; i < 18; i++) {
  278.             if (W4[i][0] == e) {
  279.                 sigma = W4[i][1];
  280.                 return sigma;
  281.             }
  282.  
  283.             if (W4[i][0] > e) {
  284.                 sigma = (W4[i - 1][1] + (W4[i][1] - W4[i - 1][1])*(e - W4[i - 1][0]) / (W4[i][0] - W4[i - 1][0]));
  285.                 return sigma;
  286.             }
  287.  
  288.         }
  289.         if (W4[17][0] < e) {
  290.             sigma = 0;
  291.             return sigma;
  292.         }
  293.         if (W4[0][0] > e) {
  294.             sigma = 0;
  295.             return sigma;
  296.         }
  297.     }
  298.     else if (type == "w5") {
  299.         for (int i = 0; i < 18; i++) {
  300.             if (W5[i][0] == e) {
  301.                 sigma = W5[i][1];
  302.                 return sigma;
  303.             }
  304.  
  305.             if (W5[i][0] > e) {
  306.                 sigma = (W5[i - 1][1] + (W5[i][1] - W5[i - 1][1])*(e - W5[i - 1][0]) / (W5[i][0] - W5[i - 1][0]));
  307.                 return sigma;
  308.             }
  309.  
  310.         }
  311.         if (W5[17][0] < e) {
  312.             sigma = 0;
  313.             return sigma;
  314.         }
  315.         if (W5[0][0] > e) {
  316.             sigma = 0;
  317.             return sigma;
  318.         }
  319.     }
  320.     else if (type == "w6") {
  321.         for (int i = 0; i < 16; i++) {
  322.             if (W6[i][0] == e) {
  323.                 sigma = W6[i][1];
  324.                 return sigma;
  325.             }
  326.  
  327.             if (W6[i][0] > e) {
  328.                 sigma = (W6[i - 1][1] + (W6[i][1] - W6[i - 1][1])*(e - W6[i - 1][0]) / (W6[i][0] - W6[i - 1][0]));
  329.                 return sigma;
  330.             }
  331.  
  332.         }
  333.         if (W6[15][0] < e) {
  334.             sigma = 0;
  335.             return sigma;
  336.         }
  337.         if (W6[0][0] > e) {
  338.             sigma = 0;
  339.             return sigma;
  340.         }
  341.     }
  342.     else if (type == "w7") {
  343.         for (int i = 0; i < 18; i++) {
  344.             if (W7[i][0] == e) {
  345.                 sigma = W7[i][1];
  346.                 return sigma;
  347.             }
  348.  
  349.             if (W7[i][0] > e) {
  350.                 sigma = (W7[i - 1][1] + (W7[i][1] - W7[i - 1][1])*(e - W7[i - 1][0]) / (W7[i][0] - W7[i - 1][0]));
  351.                 return sigma;
  352.             }
  353.  
  354.         }
  355.         if (W7[17][0] < e) {
  356.             sigma = 0;
  357.             return sigma;
  358.         }
  359.         if (W7[0][0] > e) {
  360.             sigma = 0;
  361.             return sigma;
  362.         }
  363.     }
  364.     else if (type == "w8") {
  365.         for (int i = 0; i < 16; i++) {
  366.             if (W8[i][0] == e) {
  367.                 sigma = W8[i][1];
  368.                 return sigma;
  369.             }
  370.  
  371.             if (W8[i][0] > e) {
  372.                 sigma = (W8[i - 1][1] + (W8[i][1] - W8[i - 1][1])*(e - W8[i - 1][0]) / (W8[i][0] - W8[i - 1][0]));
  373.                 return sigma;
  374.             }
  375.  
  376.         }
  377.         if (W8[15][0] < e) {
  378.             sigma = 0;
  379.             return sigma;
  380.         }
  381.         if (W8[0][0] > e) {
  382.             sigma = 0;
  383.             return sigma;
  384.         }
  385.     }
  386.     else if (type == "E1") {
  387.         for (int i = 0; i < 24; i++) {
  388.             if (CS_Ex_1[i][0] == e) {
  389.                 sigma = CS_Ex_1[i][1];
  390.                 return sigma;
  391.             }
  392.  
  393.             if (CS_Ex_1[i][0] > e) {
  394.                 sigma = (CS_Ex_1[i - 1][1] + (CS_Ex_1[i][1] - CS_Ex_1[i - 1][1])*(e - CS_Ex_1[i - 1][0]) / (CS_Ex_1[i][0] - CS_Ex_1[i - 1][0]));
  395.                 return sigma;
  396.             }
  397.  
  398.         }
  399.         if (CS_Ex_1[23][0] < e) {
  400.             sigma = 0;
  401.         }
  402.         if (CS_Ex_1[0][0] > e) {
  403.             sigma = 0;
  404.         }
  405.     }
  406.     else if (type == "E2") {
  407.         for (int i = 0; i < 24; i++) {
  408.             if (CS_Ex_2[i][0] == e) {
  409.                 sigma = CS_Ex_2[i][1];
  410.                 return sigma;
  411.             }
  412.  
  413.             if (CS_Ex_2[i][0] > e) {
  414.                 sigma = (CS_Ex_2[i - 1][1] + (CS_Ex_2[i][1] - CS_Ex_2[i - 1][1])*(e - CS_Ex_2[i - 1][0]) / (CS_Ex_2[i][0] - CS_Ex_2[i - 1][0]));
  415.                 return sigma;
  416.             }
  417.  
  418.         }
  419.         if (CS_Ex_2[23][0] < e) {
  420.             sigma = 0;
  421.         }
  422.         if (CS_Ex_2[0][0] > e) {
  423.             sigma = 0;
  424.         }
  425.     }
  426.     else if (type == "E3") {
  427.         for (int i = 0; i < 24; i++) {
  428.             if (CS_Ex_3[i][0] == e) {
  429.                 sigma = CS_Ex_3[i][1];
  430.                 return sigma;
  431.             }
  432.  
  433.             if (CS_Ex_3[i][0] > e) {
  434.                 sigma = (CS_Ex_3[i - 1][1] + (CS_Ex_3[i][1] - CS_Ex_3[i - 1][1])*(e - CS_Ex_3[i - 1][0]) / (CS_Ex_3[i][0] - CS_Ex_3[i - 1][0]));
  435.                 return sigma;
  436.             }
  437.  
  438.         }
  439.         if (CS_Ex_3[23][0] < e) {
  440.             sigma = 0;
  441.         }
  442.         if (CS_Ex_3[0][0] > e) {
  443.             sigma = 0;
  444.         }
  445.     }
  446.     else if (type == "E4") {
  447.         for (int i = 0; i < 24; i++) {
  448.             if (CS_Ex_4[i][0] == e) {
  449.                 sigma = CS_Ex_4[i][1];
  450.                 return sigma;
  451.             }
  452.  
  453.             if (CS_Ex_4[i][0] > e) {
  454.                 sigma = (CS_Ex_4[i - 1][1] + (CS_Ex_4[i][1] - CS_Ex_4[i - 1][1])*(e - CS_Ex_4[i - 1][0]) / (CS_Ex_4[i][0] - CS_Ex_4[i - 1][0]));
  455.                 return sigma;
  456.             }
  457.  
  458.         }
  459.         if (CS_Ex_4[23][0] < e) {
  460.             sigma = 0;
  461.         }
  462.         if (CS_Ex_4[0][0] > e) {
  463.             sigma = 0;
  464.         }
  465.     }
  466.     else if (type == "E5") {
  467.         for (int i = 0; i < 24; i++) {
  468.             if (CS_Ex_5[i][0] == e) {
  469.                 sigma = CS_Ex_5[i][1];
  470.                 return sigma;
  471.             }
  472.  
  473.             if (CS_Ex_5[i][0] > e) {
  474.                 sigma = (CS_Ex_5[i - 1][1] + (CS_Ex_5[i][1] - CS_Ex_5[i - 1][1])*(e - CS_Ex_5[i - 1][0]) / (CS_Ex_5[i][0] - CS_Ex_5[i - 1][0]));
  475.                 return sigma;
  476.             }
  477.  
  478.         }
  479.         if (CS_Ex_5[23][0] < e) {
  480.             sigma = 0;
  481.         }
  482.         if (CS_Ex_5[0][0] > e) {
  483.             sigma = 0;
  484.         }
  485.     }
  486.     else if (type == "E6") {
  487.         for (int i = 0; i < 24; i++) {
  488.             if (CS_Ex_6[i][0] == e) {
  489.                 sigma = CS_Ex_6[i][1];
  490.                 return sigma;
  491.             }
  492.  
  493.             if (CS_Ex_6[i][0] > e) {
  494.                 sigma = (CS_Ex_6[i - 1][1] + (CS_Ex_6[i][1] - CS_Ex_6[i - 1][1])*(e - CS_Ex_6[i - 1][0]) / (CS_Ex_6[i][0] - CS_Ex_6[i - 1][0]));
  495.                 return sigma;
  496.             }
  497.  
  498.         }
  499.         if (CS_Ex_6[23][0] < e) {
  500.             sigma = 0;
  501.         }
  502.         if (CS_Ex_6[0][0] > e) {
  503.             sigma = 0;
  504.         }
  505.     }
  506.     else if (type == "E7") {
  507.         for (int i = 0; i < 24; i++) {
  508.             if (CS_Ex_7[i][0] == e) {
  509.                 sigma = CS_Ex_7[i][1];
  510.                 return sigma;
  511.             }
  512.  
  513.             if (CS_Ex_7[i][0] > e) {
  514.                 sigma = (CS_Ex_7[i - 1][1] + (CS_Ex_7[i][1] - CS_Ex_7[i - 1][1])*(e - CS_Ex_7[i - 1][0]) / (CS_Ex_7[i][0] - CS_Ex_7[i - 1][0]));
  515.                 return sigma;
  516.             }
  517.  
  518.         }
  519.         if (CS_Ex_7[23][0] < e) {
  520.             sigma = 0;
  521.         }
  522.         if (CS_Ex_7[0][0] > e) {
  523.             sigma = 0;
  524.         }
  525.     }
  526.     else if (type == "E8") {
  527.         for (int i = 0; i < 18; i++) {
  528.             if (CS_Ex_8[i][0] == e) {
  529.                 sigma = CS_Ex_8[i][1];
  530.                 return sigma;
  531.             }
  532.  
  533.             if (CS_Ex_8[i][0] > e) {
  534.                 sigma = (CS_Ex_8[i - 1][1] + (CS_Ex_8[i][1] - CS_Ex_8[i - 1][1])*(e - CS_Ex_8[i - 1][0]) / (CS_Ex_8[i][0] - CS_Ex_8[i - 1][0]));
  535.                 return sigma;
  536.             }
  537.  
  538.         }
  539.         if (CS_Ex_8[17][0] < e) {
  540.             sigma = 0;
  541.         }
  542.         if (CS_Ex_8[0][0] > e) {
  543.             sigma = 0;
  544.         }
  545.     }
  546.     else if (type == "E9") {
  547.         for (int i = 0; i < 28; i++) {
  548.             if (CS_Ex_9[i][0] == e) {
  549.                 sigma = CS_Ex_9[i][1];
  550.                 return sigma;
  551.             }
  552.  
  553.             if (CS_Ex_9[i][0] > e) {
  554.                 sigma = (CS_Ex_9[i - 1][1] + (CS_Ex_9[i][1] - CS_Ex_9[i - 1][1])*(e - CS_Ex_9[i - 1][0]) / (CS_Ex_9[i][0] - CS_Ex_9[i - 1][0]));
  555.                 return sigma;
  556.             }
  557.  
  558.         }
  559.         if (CS_Ex_9[27][0] < e) {
  560.             sigma = 0;
  561.         }
  562.         if (CS_Ex_9[0][0] > e) {
  563.             sigma = 0;
  564.         }
  565.     }
  566.     else if (type == "E10") {
  567.         for (int i = 0; i < 16; i++) {
  568.             if (CS_Ex_10[i][0] == e) {
  569.                 sigma = CS_Ex_10[i][1];
  570.                 return sigma;
  571.             }
  572.  
  573.             if (CS_Ex_10[i][0] > e) {
  574.                 sigma = (CS_Ex_10[i - 1][1] + (CS_Ex_10[i][1] - CS_Ex_10[i - 1][1])*(e - CS_Ex_10[i - 1][0]) / (CS_Ex_10[i][0] - CS_Ex_10[i - 1][0]));
  575.                 return sigma;
  576.             }
  577.  
  578.         }
  579.         if (CS_Ex_10[15][0] < e) {
  580.             sigma = 0;
  581.         }
  582.         if (CS_Ex_10[0][0] > e) {
  583.             sigma = 0;
  584.         }
  585.     }
  586.     else if (type == "E11") {
  587.         for (int i = 0; i < 24; i++) {
  588.             if (CS_Ex_11[i][0] == e) {
  589.                 sigma = CS_Ex_11[i][1];
  590.                 return sigma;
  591.             }
  592.  
  593.             if (CS_Ex_11[i][0] > e) {
  594.                 sigma = (CS_Ex_11[i - 1][1] + (CS_Ex_11[i][1] - CS_Ex_11[i - 1][1])*(e - CS_Ex_11[i - 1][0]) / (CS_Ex_11[i][0] - CS_Ex_11[i - 1][0]));
  595.                 return sigma;
  596.             }
  597.  
  598.         }
  599.         if (CS_Ex_11[23][0] < e) {
  600.             sigma = 0;
  601.         }
  602.         if (CS_Ex_11[0][0] > e) {
  603.             sigma = 0;
  604.         }
  605.     }
  606.  
  607.  
  608.     return sigma;
  609. }
  610.  
  611. void main() {
  612.  
  613.  
  614.  
  615.     double Cmaxwell = pow((2*PI*m*T0),-1.5);   // константа в функции распределения максвелла
  616.  
  617.  
  618.     init();
  619.  
  620.     double **C_S = new double*[13];
  621.  
  622.     for (int i = 0; i < 21; i++) {
  623.         C_S[i] = new double[M];
  624.     }
  625.  
  626.     for (int j = 0; j < M; j++) {
  627.         C_S[0][j] = cs(h*j, "t");
  628.         C_S[1][j] = cs(h*j, "r");
  629.         C_S[2][j] = cs(h*j, "E1");
  630.         C_S[3][j] = cs(h*j, "E2");
  631.         C_S[4][j] = cs(h*j, "E3");
  632.         C_S[5][j] = cs(h*j, "E4");
  633.         C_S[6][j] = cs(h*j, "E5");
  634.         C_S[7][j] = cs(h*j, "E6");
  635.         C_S[8][j] = cs(h*j, "E7");
  636.         C_S[9][j] = cs(h*j, "E8");
  637.         C_S[10][j] = cs(h*j, "E9");
  638.         C_S[11][j] = cs(h*j, "E10");
  639.         C_S[12][j] = cs(h*j, "E11");
  640.         C_S[13][j] = cs(h*j, "w1");
  641.         C_S[14][j] = cs(h*j, "w2");
  642.         C_S[15][j] = cs(h*j, "w3");
  643.         C_S[16][j] = cs(h*j, "w4");
  644.         C_S[17][j] = cs(h*j, "w5");
  645.         C_S[18][j] = cs(h*j, "w6");
  646.         C_S[19][j] = cs(h*j, "w7");
  647.         C_S[20][j] = cs(h*j, "w8");
  648.     }
  649.  
  650.  
  651.     FILE *FREE; FREE = fopen("C:\\FREE.txt", "w");
  652.     FILE *integral; integral = fopen("C:\\integral.txt", "w");
  653.  
  654.     double A, B, C, D, E, Lambda;                                                       // константы в изначальной разностной схеме
  655.     double C1, C2, C3, C4;                                                      // контанты в конечной разностной схеме
  656.  
  657.     double I = 0;                                                // нормировочный интеграл
  658.  
  659.     double Fp[M], Fn[M]; // Fp - предыдуший шаг по времени, Fn - следующий шаг по времени
  660.  
  661.     double F[z][M];             // массив куда записывается функция распределения на нужных шагах
  662.  
  663.  
  664. //  пороги возбуждения электронных уровней
  665.     int Ts_Ex_1 = (int)floor(6.17 / h);
  666.     int Ts_Ex_2 = (int)floor(7 / h);
  667.     int Ts_Ex_3 = (int)floor(7.35 / h);
  668.     int Ts_Ex_4 = (int)floor(7.36 / h);
  669.     int Ts_Ex_5 = (int)floor(8.16 / h);
  670.     int Ts_Ex_6 = (int)floor(8.4 / h);
  671.     int Ts_Ex_7 = (int)floor(8.55 / h);
  672.     int Ts_Ex_8 = (int)floor(8.89 / h);
  673.     int Ts_Ex_9 = (int)floor(11.03 / h);
  674.     int Ts_Ex_10 = (int)floor(11.87 / h);
  675.     int Ts_Ex_11 = (int)floor(12.25 / h);
  676.  
  677.     int TS_Vibr_01 = (int)floor(0.3 / h);
  678.     int TS_Vibr_02 = (int)floor(1.8 / h);
  679.     int TS_Vibr_03 = (int)floor(2 / h);
  680.     int TS_Vibr_04 = (int)floor(2.1 / h);
  681.     int TS_Vibr_05 = (int)floor(2.2 / h);
  682.     int TS_Vibr_06 = (int)floor(2.3 / h);
  683.     int TS_Vibr_07 = (int)floor(2.4 / h);
  684.     int TS_Vibr_08 = (int)floor(2.6 / h);
  685.  
  686.     double SUM_CS_VIBR = 0;               //сумма сечений колебательного возбуждения
  687.     double VIBR_1_01, VIBR_1_02, VIBR_1_03, VIBR_1_04, VIBR_1_05, VIBR_1_06, VIBR_1_07, VIBR_1_08 ;//константа для прямых ударов второго рода
  688.     double VIBR_2_01, VIBR_2_02, VIBR_2_03, VIBR_2_04, VIBR_2_05, VIBR_2_06, VIBR_2_07, VIBR_2_08 ; //константы для обратных ударов второго рода
  689.     double SUM_VIBR1 = 0;   // сумма констант для прямых ударов второго рода
  690.     double SUM_VIBR2 = 0; // сумма констант для обратных ударов первого рода
  691.     double SUM_VIBR3 = 0; //сумма для констант обратных ударов второго рода
  692.  
  693.     double G, G1, G2, G3, G4, G5, G6, G7, G8, G9, G10, G11;  //отн. засел. колеб. сост
  694.     double G0; //отн. засел. эл. сост
  695.  
  696.     double Ex_sum = 0; //сумма сечений колебательного возбуждения
  697.     double Ex1 = 0, Ex2 = 0, Ex3 = 0, Ex4 = 0, Ex5 = 0, Ex6 = 0, Ex7 = 0, Ex8 = 0, Ex9 = 0, Ex10 = 0, Ex11 = 0;     //константа колебательного возбуждения для ударов второго рода
  698.     double SUM_EX = 0; // сумма констант прямых ударов второго рода
  699.  
  700.     double PartFun = 0.847037;
  701.  
  702.    
  703.  
  704.     A = sqrt(m / 2)*sqrt(eps0)/n0;
  705.     B = ((q*En*1e-19)*(q*En*1e-19))/ (9*3 * eps0);
  706.     C = d*T0*eps0;
  707.     D = d*eps0;
  708.     E = 4 * B0*eps0;
  709.  
  710.     G = eps0*exp((hw0*pow((0 + 0.5),2) - hw*(0 + 0.5))/Tk)/ PartFun;
  711.     G1 = eps0*exp((hw0*pow((1 + 0.5), 2) - hw*(1 + 0.5)) / Tk) / PartFun;
  712.     G2 = eps0*exp((hw0*pow((2 + 0.5), 2) - hw*(2 + 0.5)) / Tk) / PartFun;
  713.     G3 = eps0*exp((hw0*pow((3 + 0.5), 2) - hw*(3 + 0.5)) / Tk) / PartFun;
  714.     G4 = eps0*exp((hw0*pow((4 + 0.5), 2) - hw*(4 + 0.5)) / Tk) / PartFun;
  715.     G5 = eps0*exp((hw0*pow((5 + 0.5), 2) - hw*(5 + 0.5)) / Tk) / PartFun;
  716.     G6 = eps0*exp((hw0*pow((6 + 0.5), 2) - hw*(6 + 0.5)) / Tk) / PartFun;
  717.     G7 = eps0*exp((hw0*pow((7 + 0.5), 2) - hw*(7 + 0.5)) / Tk) / PartFun;
  718.     G8 = eps0*exp((hw0*pow((8 + 0.5), 2) - hw*(8 + 0.5)) / Tk) / PartFun;
  719.     G9 = eps0*exp((hw0*pow((9 + 0.5), 2) - hw*(9 + 0.5)) / Tk) / PartFun;
  720.     G10 = eps0*exp((hw0*pow((10 + 0.5), 2) - hw*(10 + 0.5)) / Tk) / PartFun;
  721.     G11 = eps0*exp((hw0*pow((11 + 0.5), 2) - hw*(11 + 0.5)) / Tk) / PartFun;
  722.     G0 = eps0 * 1 / 1.07674;   //Электронные состояния
  723.  
  724.  
  725.     Lambda = 2 * PI*pow(q, 4)*(ne / n0)*logL;  
  726.  
  727.    
  728.     for (int j = 0; j < M; j++) {   //начальное условие
  729.         Fp[j] = 2 * PI*pow(2 * m, 1.5)*Cmaxwell*exp(-h*j / T0);
  730.  
  731.  
  732.     }
  733.  
  734.  
  735.  
  736.     for (int i = 1; i < N; i++) {
  737.         for (int j = 1; j < M - 1; j++) {
  738.  
  739.  
  740.  
  741.             C1 = (B*(h*j) / C_S[0][j]) + C*(h*j)*(h*j)*C_S[0][j];
  742.             C2 = D*(h*j)*(h*j)*C_S[0][j] + E*(h*j)*C_S[1][j];
  743.             C3 = B*(C_S[0][j] - j*(C_S[0][j+1] - C_S[0][j])) / (C_S[0][j] * C_S[0][j]) +2 * C*(h*j)*C_S[0][j] + C*(h*j)*(j)*(C_S[0][j+1] - C_S[0][j]) + C2;
  744.             C4 = 2 * D*(h*j)*C_S[0][j] + D*(h*j)*(j)*(C_S[0][j+1] - C_S[0][j]) + E*C_S[1][j] + E*(j)*(C_S[1][j+1] - C_S[1][j]);
  745.  
  746.             if (j + Ts_Ex_1 < M) {
  747.                 Ex1 = (h*j + 6.17)*C_S[2][j + Ts_Ex_1] * Fp[j + Ts_Ex_1];
  748.             }
  749.             else (Ex1 = 0);
  750.  
  751.             if (j + Ts_Ex_2 < M) {
  752.                 Ex2 = (h*j + 7)*C_S[3][j + Ts_Ex_2] * Fp[j + Ts_Ex_2];
  753.             }
  754.             else (Ex2 = 0);
  755.  
  756.             if (j + Ts_Ex_3 < M) {
  757.                 Ex3 = (h*j + 7.35)*C_S[4][j + Ts_Ex_3] * Fp[j + Ts_Ex_3];
  758.             }
  759.             else (Ex3 = 0);
  760.  
  761.             if (j + Ts_Ex_4 < M) {
  762.                 Ex4 = (h*j + 7.36)*C_S[5][j + Ts_Ex_4] * Fp[j + Ts_Ex_4];
  763.             }
  764.             else (Ex4 = 0);
  765.  
  766.             if (j + Ts_Ex_5 < M) {
  767.                 Ex5 = (h*j + 8.16)*C_S[6][j + Ts_Ex_5] * Fp[j + Ts_Ex_5];
  768.             }
  769.             else (Ex5 = 0);
  770.  
  771.             if (j + Ts_Ex_6 < M) {
  772.                 Ex6 = (h*j + 8.4)*C_S[7][j + Ts_Ex_6] * Fp[j + Ts_Ex_6];
  773.             }
  774.             else (Ex6 = 0);
  775.  
  776.             if (j + Ts_Ex_7 < M) {
  777.                 Ex7 = (h*j + 8.55)*C_S[8][j + Ts_Ex_7] * Fp[j + Ts_Ex_7];
  778.             }
  779.             else (Ex7 = 0);
  780.  
  781.             if (j + Ts_Ex_8 < M) {
  782.                 Ex8 = (h*j + 8.89)*C_S[9][j + Ts_Ex_8] * Fp[j + Ts_Ex_8];
  783.             }
  784.             else (Ex8 = 0);
  785.  
  786.             if (j + Ts_Ex_9 < M) {
  787.                 Ex9 = (h*j + 11.03)*C_S[10][j + Ts_Ex_9] * Fp[j + Ts_Ex_9];
  788.             }
  789.             else (Ex9 = 0);
  790.  
  791.             if (j + Ts_Ex_10 < M) {
  792.                 Ex10 = (h*j + 11.87)*C_S[11][j + Ts_Ex_10] * Fp[j + Ts_Ex_10];
  793.             }
  794.             else (Ex10 = 0);
  795.  
  796.             if (j + Ts_Ex_11 < M) {
  797.                 Ex11 = (h*j + 12.25)*C_S[12][j + Ts_Ex_11] * Fp[j + Ts_Ex_11];
  798.             }
  799.             else (Ex11 = 0);
  800.  
  801.             SUM_EX = Ex1 + Ex2 + Ex3 + Ex4 + Ex5 + Ex6 + Ex7 + Ex8 + Ex9 + Ex10 + Ex11;
  802.  
  803.  
  804.             Ex_sum = C_S[2][j] + C_S[3][j] + C_S[4][j] + C_S[5][j] + C_S[6][j] + C_S[7][j] + C_S[8][j] + C_S[9][j] + C_S[10][j] + C_S[11][j] + C_S[12][j];
  805.  
  806.  
  807.             if (j + TS_Vibr_01 < M) {
  808.                 VIBR_1_01 = (h*j + 0.3)*C_S[13][j + TS_Vibr_01];
  809.             }
  810.             else (VIBR_1_01 = 0);
  811.  
  812.             if (j + TS_Vibr_02 < M) {
  813.                 VIBR_1_02 = (h*j + 1.8)*C_S[14][j + TS_Vibr_02];
  814.             }
  815.             else (VIBR_1_02 = 0);
  816.  
  817.             if (j + TS_Vibr_03 < M) {
  818.                 VIBR_1_03 = (h*j + 2)*C_S[15][j + TS_Vibr_03];
  819.             }
  820.             else (VIBR_1_03 = 0);
  821.  
  822.             if (j + TS_Vibr_04 < M) {
  823.                 VIBR_1_04 = (h*j + 2.1)*C_S[16][j + TS_Vibr_04];
  824.             }
  825.             else (VIBR_1_04 = 0);
  826.  
  827.             if (j + TS_Vibr_05 < M) {
  828.                 VIBR_1_05 = (h*j + 2.2)*C_S[17][j + TS_Vibr_05];
  829.             }
  830.             else (VIBR_1_05 = 0);
  831.  
  832.             if (j + TS_Vibr_06 < M) {
  833.                 VIBR_1_06 = (h*j + 2.3)*C_S[18][j + TS_Vibr_06];
  834.             }
  835.             else (VIBR_1_06 = 0);
  836.  
  837.             if (j + TS_Vibr_07 < M) {
  838.                 VIBR_1_07 = (h*j + 2.4)*C_S[19][j + TS_Vibr_07];
  839.             }
  840.             else (VIBR_1_07 = 0);
  841.  
  842.             if (j + TS_Vibr_08 < M) {
  843.                 VIBR_1_08 = (h*j + 2.6)*C_S[20][j + TS_Vibr_08];
  844.             }
  845.             else (VIBR_1_08 = 0);
  846.  
  847.             if (j - TS_Vibr_01 >= 0) {
  848.                 VIBR_2_01 = C_S[13][j] * Fp[j - TS_Vibr_01];
  849.             }
  850.             else (VIBR_2_01 = 0);
  851.  
  852.             if (j - TS_Vibr_02 >= 0) {
  853.                 VIBR_2_02 = C_S[14][j] * Fp[j - TS_Vibr_02];
  854.             }
  855.             else (VIBR_2_02 = 0);
  856.  
  857.             if (j - TS_Vibr_03 >= 0) {
  858.                 VIBR_2_03 = C_S[15][j] * Fp[j - TS_Vibr_03];
  859.             }
  860.             else (VIBR_2_03 = 0);
  861.  
  862.             if (j - TS_Vibr_04 >= 0) {
  863.                 VIBR_2_04 = C_S[16][j] * Fp[j - TS_Vibr_04];
  864.             }
  865.             else (VIBR_2_04 = 0);
  866.  
  867.             if (j - TS_Vibr_05 >= 0) {
  868.                 VIBR_2_05 = C_S[17][j] * Fp[j - TS_Vibr_05];
  869.             }
  870.             else (VIBR_2_05 = 0);
  871.  
  872.             if (j - TS_Vibr_06 >= 0) {
  873.                 VIBR_2_06 = C_S[18][j] * Fp[j - TS_Vibr_06];
  874.             }
  875.             else (VIBR_2_06 = 0);
  876.  
  877.             if (j - TS_Vibr_07 >= 0) {
  878.                 VIBR_2_07 = C_S[19][j] * Fp[j - TS_Vibr_07];
  879.             }
  880.             else (VIBR_2_07 = 0);
  881.  
  882.             if (j - TS_Vibr_08 >= 0) {
  883.                 VIBR_2_08 = C_S[20][j] * Fp[j - TS_Vibr_08];
  884.             }
  885.             else (VIBR_2_08 = 0);
  886.  
  887.  
  888.  
  889.             SUM_VIBR1 = VIBR_1_01 * Fp[j + TS_Vibr_01] + VIBR_1_02 * Fp[j + TS_Vibr_02] + VIBR_1_03 * Fp[j + TS_Vibr_03] +
  890.                 VIBR_1_04 * Fp[j + TS_Vibr_04] + VIBR_1_05 * Fp[j + TS_Vibr_05] + VIBR_1_06 * Fp[j + TS_Vibr_06] + VIBR_1_07 * Fp[j + TS_Vibr_07] + VIBR_1_08 * Fp[j + TS_Vibr_08];
  891.  
  892.             SUM_CS_VIBR = C_S[13][j] + C_S[14][j] + C_S[15][j] + C_S[16][j] + C_S[17][j] + C_S[18][j] + C_S[19][j] + C_S[20][j];
  893.  
  894.             SUM_VIBR2 = G1*VIBR_1_01*exp(hw/Tk) + G2*VIBR_1_02*exp(2*hw/Tk) + G3*VIBR_1_03*exp(3*hw/Tk) + G4*VIBR_1_04*exp(4*hw/Tk) + G5*VIBR_1_05*exp(5*hw/Tk) + G6*VIBR_1_06*exp(6*hw/Tk) + G7*VIBR_1_07*exp(7*hw/Tk) + G8*VIBR_1_08*exp(8*hw/Tk);
  895.  
  896.             SUM_VIBR3 = G1*VIBR_2_01*exp(hw/Tk) + G2*VIBR_2_02*exp(2*hw/Tk) + G3*VIBR_2_03*exp(3*hw/Tk) + G4*VIBR_2_04*exp(4*hw/Tk) + G5*VIBR_2_05*(5*hw/Tk) + G6*VIBR_2_06*(6*hw/Tk) + G7*VIBR_2_07*exp(7*hw/Tk) + G8*VIBR_2_08*exp(8*hw/Tk);
  897.                        
  898.  
  899.             Fn[j] = Fp[j] + (t / (A*sqrt(h*j))) * ((Fp[j+1]+Fp[j-1] - 2*Fp[j])*C1/(h*h) + (Fp[j]-Fp[j-1])*C3/h + Fp[j]*C4
  900.                 - G0*(h*j)*Fp[j] * Ex_sum + G0*SUM_EX
  901.                 - G*(h*j)*Fp[j] * SUM_CS_VIBR + G*SUM_VIBR1
  902.                 -SUM_VIBR2*Fp[j] + (h*j)*SUM_VIBR3);
  903.  
  904.            
  905.             if (i == 2) { F[0][j] = Fn[j]; }
  906.             if (i == 10) { F[1][j] = Fn[j]; }
  907.             if (i == 100) { F[2][j] = Fn[j]; }
  908.             if (i == 1000) { F[3][j] = Fn[j]; }
  909.             if (i == 10000) { F[4][j] = Fn[j]; }
  910.             if (i == 100000) { F[5][j] = Fn[j]; }
  911.             if (i == 1000000) { F[6][j] = Fn[j]; }
  912.             if (i == 10000000) { F[7][j] = Fn[j]; }
  913. //          if (i == 100000000) { F[8][j] = Fn[j]; }
  914.    
  915.             I = I + Fn[j] * sqrt(h*j)*h;
  916.  
  917.  
  918.  
  919. //          printf("%e\n", SUM_VIBR3);
  920. //          Sleep(0);
  921.         }
  922. //      printf("\n");
  923. //      Sleep(1000);
  924.        
  925.         Fn[M - 1] = 0;      // Добавили граничные условия справа
  926.         Fn[0] = Fn[1];      // и слева
  927.        
  928.  
  929. //      I = I * (2 * PI)*pow((2 * m), 1.5);
  930.  
  931.  
  932.             for (int j = 0; j < M; j++) {
  933.                 Fp[j] = Fn[j];
  934.             }
  935. //          printf("%e\n", I);
  936. //          Sleep(10000);
  937.  
  938.             I = 0;
  939.  
  940.        
  941.     }
  942.  
  943.  
  944.     for (int j = 0; j < M; j++) {
  945.         for (int i = 0; i < z; i = i++) {
  946.             F[i][0] = F[i][1];
  947.             F[i][M - 1] = 0;
  948.  
  949.             fprintf(FREE, "%e\t", F[i][j]);
  950.         }
  951.         fprintf(FREE, "\n");
  952.     }
  953.    
  954.     for (int i = 0; i < 2; i++) {
  955.         delete[]C_S[i];
  956.     }
  957.     fclose(FREE);
  958.     fclose(integral);
  959.  
  960.  
  961.        
  962. }
Add Comment
Please, Sign In to add comment
Public Pastes
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy.  OK, I Understand
Not a member of Pastebin yet?
Sign Up, it unlocks many cool features!