Редактирование: Моделирование углекислого газа методом динамики частиц

[[ Курсовые_работы_по_ВМДС:_2016-2017 | Курсовые работы 2016-2017 учебного года]] > '''Моделирование углекислого газа методом динамики частиц''' <HR>

'''''Курсовой проект по [[Механика дискретных сред|Механике дискретных сред]]'''''

'''Исполнитель:''' [[Смирнов Александр]]

'''Группа:''' [[Группа 09|09]] (43604/1)

'''Семестр:''' осень 2016


== Формулировка задачи ==
[[Файл:foto11.jpg|thumb|Рис.1 Молекула углекислого газа.|400px]]
Смоделировать молекулы углекислого газа методом динамики частиц, проверить выполнение закона сохранения энергии и рассмотреть распределение энергии по степеням свободы.

== Общие сведения ==
Для моделирования используем частицы, которые представляют собой абсолютно упругий шар. Масса углерода равна 12, а кислорода - 16 условных единиц. Для отталкивания молекул используется потенциал Морзе, внутри молекулы: Упругая сила и угловая пружина.
===== Упругая сила и угловая пружина =====
Упругая сила определяется формулой:
::<math>
    F( r) = -kr,
</math>
где
* <math>k</math> — жесткость,
* <math> r</math> — отклонение от положения равновесия.
Потенциал пружины равен:
::<math>
    \varPi(r) =  С\frac{(\pi-φ)^2}{2}
    ,
</math>
При этом
::<math>φ =\frac{ {\bf r_1}·{\bf r_2}}{r_1 r_2} </math>
где
* <math>С</math> —угловая жесткость,
* <math> {\bf r_i} </math> — радиус вектор от атома углерода к первому или второму кислороду соответственно.
А сила вычисляется как:
::<math>F(r) = -\varPi'(r) ,
</math>
 

===== [[Потенциал Морзе |  Потенциал Морзе ]] =====
Парный силовой потенциал взаимодействия.
Определяется формулой: 
::<math>
    \varPi(r) = 
    D\left[e^{-2\alpha(r-a)}-2e^{-\alpha(r-a)}\right],
</math>

где
* <math>D</math> — энергия связи,
* <math>a</math> — длина связи,
* <math>\alpha</math> — параметр, характеризующий ширину потенциальной ямы.

Сила, соответствующая потенциалу Морзе, вычисляется по формуле:
::<math>
    F(r) = -\varPi'(r) = 
    2\alpha D\left[e^{-2\alpha(r-a)}-e^{-\alpha(r-a)}\right].
</math>
Или в векторной форме: 
::<math>
    {\bf F}({\bf r})= -\nabla\varPi(r) = 2\alpha D\left[e^{-2\alpha(r-a)}-e^{-\alpha(r-a)}\right]\frac{{\bf r}}{r}
</math>
==Программа==
В данной программе в начальный момент времени системе задаются случайные скорости(начальная энергия,они достаточно велики, чтобы можно было пренебречь потенциальной энергией взаимодействия. Можно менять количество молекул углекислого газа и сбрасывать таймер расчета средних значений. Так же выводятся:кинетическая энергия системы в данный момент времени, средняя кинетическая энергия системы в данный момент времени, средние энергии, приходящиеся на атом углерода, первый и второй атом кислорода в молекуле. 
{{#widget:Iframe |url=https://ailurus.ru/protos/foten/lab-k/ |width=680 |height=680 |border=0 }} 

<div class="mw-collapsible mw-collapsed" style="width:100%" >
'''Текст программы на языке JavaScript:'''
<div class="mw-collapsible-content">
Файл '''"Lab.js"'''
<syntaxhighlight lang="javascript" line start="1" enclose="div">
window.addEventListener("load", MainBalls, true);
function MainBalls() {

    // Предварительные установки

    var canvas = canvasBalls;
    var context = canvas.getContext("2d");                  // на context происходит рисование
    canvas.oncontextmenu = function (e) {return false;};    // блокировка контекстного меню

    var Pi = 3.1415926;                   // число "пи"

    var m0 = 1;                           // масштаб массы
    var T0 = 1;                           // масштаб времени (период колебаний исходной системы)
    var a0 = 1;                           // масштаб расстояния (диаметр шара)

    var g0 = a0 / T0 / T0;                // масштаб ускорения (ускорение, при котором за T0 будет пройдено расстояние a0)
    var k0 = 2 * Pi / T0;                 // масштаб частоты
    var C0 = m0 * k0 * k0;                // масштаб жесткости
    var B0 = 2 * m0 * k0;                 // масштаб вязкости

    // *** Задание физических параметров ***

    var Ny = 5;                           // число шаров, помещающихся по вертикали в окно (задает размер шара относительно размера окна)
    var m = 12 * m0;                       // масса C
	var m1 = 16 * m0 ;                   //масса О 
    var Cwall = 8 * C0;                  // жесткость стен
    var B = 0 * B0;                   // вязкость среды
    var Bwall = 0 * B0;                // вязкость на стенках
    var Cball =  Cwall;              // жесткость между частицами            
    var r = 0.1 * a0;                     // радиус частицы в расчетных координатах
	var K = 0.85;                         // сила взаимодействия ограничивается значением, реализующимся при r/a = K
    var a = 2.2 * r;                        // равновесное расстояние между частицами
	var a_1 = 2.05 * r;                       
    var a_2 = 2.1 * a ;
	var aCut =  2.05 * r ;                     // радиус обрезания
	var alfa = 4 ;
	var j = 20 ;
	var Ek = 0 ;
	var Ek1 = 0 ;
	var k = 1 ;
	var Ekc = 0 ;
	var Ekck = 0;
	var EC = 0 ;
	var EO1 = 0;
	var EO2 = 0;
	var ECk = 0 ;
	var EO1k = 0;
	var EO2k = 0;

    // *** Задание вычислительных параметров ***

    var fps = 60;                         // frames per second - число кадров в секунду (качечтво отображения)
    var spf = 400;                        // steps per frame   - число шагов интегрирования между кадрами (скорость расчета)
    var dt  = 0.03 * T0 / fps;           // шаг интегрирования (качество расчета)

    // Выполнение программы

    var scale = canvas.height / Ny / a0 ;  // масштабный коэффициент для перехода от расчетных к экранным координатам
    var r2 = r * r ;                       // ___в целях оптимизации___
    var aCut2 = aCut * aCut ;              // ___в целях оптимизации___
    var a2 = a * a ;                       // ___в целях оптимизации___
    var a22 = a_2 * a_2 ;
    var a11 = a_1 * a_1 ;  	
	var D = a2 * Cball / 72 ;              // энергия связи между частицами
    var LJCoeff = 12 * D / a2 ;            // коэффициент для расчета потенциала Л-Дж
	var a1 = alfa / a ;
	var MorzCoeff = 2 * a1 * D  ; 
 
    var Ka = K * r ;                       // ___в целях оптимизации___
    var K2a2 = Ka*Ka ;                // ___в целях оптимизации___

    var w = canvas.width / scale ;           // ширина окна в расчетных координатах
    var h = canvas.height / scale ;          // высота окна в расчетных координатах
	
	
	 var dNd = null ;                         // ссылка на захваченный курсором шар (drag & drop)

    // Работа с мышью

    var mx_, my_;                           // буфер позиции мыши (для расчета скорости при отпускании шара)

    canvas.onmousedown = function(e) {      // функция при нажатии клавиши мыши
        var m = mouseCoords(e);             // получаем расчетные координаты курсора мыши
        // цикл в обратную сторону, чтобы захватывать шар, нарисованный "сверху"
        // (т.к. цикл рисования идет в обычном порядке)
        for (var i = balls.length - 1; i >= 0; i--) {
            var b = balls[i];
            var rx = b.x - m.x;
            var ry = b.y - m.y;
            var rLen2 = rx * rx + ry * ry;      // квадрат расстояния между курсором и центром шара
            if (rLen2 <= r2) {                  // курсор нажал на шар
                if (e.which == 1) {             // нажата левая клавиша мыши
                    dNd = b;
                    dNd.xPlus = dNd.x - m.x;    // сдвиг курсора относительно центра шара по x
                    dNd.yPlus = dNd.y - m.y;    // сдвиг курсора относительно центра шара по y
                    mx_ = m.x;    my_ = m.y;
                    canvas.onmousemove = mouseMove;     // пока клавиша нажата - работает функция перемещения
                } else if (e.which == 3) {              // нажата правая клавиша мыши
                    balls.splice(i, 1);                 // удалить шар
                }
                return;
            }
        }
		for (var i = C.length - 1; i >= 0; i--) {
            var b = C[i];
            var rx = b.x - m.x;
            var ry = b.y - m.y;
            var rLen2 = rx * rx + ry * ry;      // квадрат расстояния между курсором и центром шара
            if (rLen2 <= r2) {                  // курсор нажал на шар
                if (e.which == 1) {             // нажата левая клавиша мыши
                    dNd = b;
                    dNd.xPlus = dNd.x - m.x;    // сдвиг курсора относительно центра шара по x
                    dNd.yPlus = dNd.y - m.y;    // сдвиг курсора относительно центра шара по y
                    mx_ = m.x;    my_ = m.y;
                    canvas.onmousemove = mouseMove;     // пока клавиша нажата - работает функция перемещения
                } else if (e.which == 3) {              // нажата правая клавиша мыши
                    C.splice(i, 1);                 // удалить шар
                }
                return;
            }
        }
		for (var i = O1.length - 1; i >= 0; i--) {
            var b = O1[i];
            var rx = b.x - m.x;
            var ry = b.y - m.y;
            var rLen2 = rx * rx + ry * ry;      // квадрат расстояния между курсором и центром шара
            if (rLen2 <= r2) {                  // курсор нажал на шар
                if (e.which == 1) {             // нажата левая клавиша мыши
                    dNd = b;
                    dNd.xPlus = dNd.x - m.x;    // сдвиг курсора относительно центра шара по x
                    dNd.yPlus = dNd.y - m.y;    // сдвиг курсора относительно центра шара по y
                    mx_ = m.x;    my_ = m.y;
                    canvas.onmousemove = mouseMove;     // пока клавиша нажата - работает функция перемещения
                } else if (e.which == 3) {              // нажата правая клавиша мыши
                    O1.splice(i, 1);                 // удалить шар
                }
                return;
            }
        }
		for (var i = O2.length - 1; i >= 0; i--) {
            var b = O2[i];
            var rx = b.x - m.x;
            var ry = b.y - m.y;
            var rLen2 = rx * rx + ry * ry;      // квадрат расстояния между курсором и центром шара
            if (rLen2 <= r2) {                  // курсор нажал на шар
                if (e.which == 1) {             // нажата левая клавиша мыши
                    dNd = b;
                    dNd.xPlus = dNd.x - m.x;    // сдвиг курсора относительно центра шара по x
                    dNd.yPlus = dNd.y - m.y;    // сдвиг курсора относительно центра шара по y
                    mx_ = m.x;    my_ = m.y;
                    canvas.onmousemove = mouseMove;     // пока клавиша нажата - работает функция перемещения
                } else if (e.which == 3) {              // нажата правая клавиша мыши
                    O2.splice(i, 1);                 // удалить шар
                }
                return;
            }
        }

        // если не вышли по return из цикла - нажатие было вне шара, добавляем новый
       if (e.which == 1) {
            dNd = addNewBall(m.x, m.y);         // добавляем шар и сразу захватываем его курсором
            if (dNd == null) return;            // если шар не добавился (из за стен или других шаров) - возвращаемся
            dNd.xPlus = 0;  dNd.yPlus = 0;      // держим шар по центру
            mx_ = m.x;    my_ = m.y;
            canvas.onmousemove = mouseMove;     // пока клавиша нажата - работает функция перемещения
        }

	
	};
	this.setNy = function(ny) {
        j = ny;   
    };
	this.ResetTime = function() {
        timers();   
    };
	
	this.newSystem = function() {
		 balls = [];  
		C =[] ; 
		O1 = [] ;
		O2 = [] ;
		if(j == 2){
		r = 0.2 * a0 ;
		 a = 3.1 * r;                        // равновесное расстояние между частицами
	 a_1 = 2.1 * r;                       
     a_2 = 2.1 * a ;
	 aCut =  2 * r ; 
	 
	
		
		r2 = r * r ;                      
		aCut2 = aCut * aCut ;              
		a2 = a * a ;                       
		a22 = a_2 * a_2 ;
		a11 = a_1 * a_1 ;  	
		D = a2 * Cball / 72 ;              
		LJCoeff = 12 * D / a2 ;            
		a1 = alfa / a ;
		MorzCoeff = 2 * a1 * D  ; 
		Ka = K * r ;                       
		K2a2 = Ka*Ka ;        
		
		          
       
	   
				addNewC( w* 0.25 , h*0.2, 10 );
				addNewO1(w* 0.25 + a , h*0.2 , 10 );
				addNewO2(w* 0.25 - a, h*0.2 , 10 );  
		  
				addNewC( w* 0.75 , h*0.2, 10 );
				addNewO1(w* 0.75 + a ,h*0.2, 10 );
				addNewO2(w* 0.75 - a, h*0.2, 10 );  
		   
		
	}
		
		if(j == 8){
		r = 0.2 * a0 ;
		 a = 2.5 * r;                        // равновесное расстояние между частицами
	 a_1 = 2.1 * r;                       
     a_2 = 2.1 * a ;
	 aCut =  2 * r ; 
	 
		j = j*0.5 ;
		
		r2 = r * r ;                      
		aCut2 = aCut * aCut ;              
		a2 = a * a ;                       
		a22 = a_2 * a_2 ;
		a11 = a_1 * a_1 ;  	
		D = a2 * Cball / 72 ;              
		LJCoeff = 12 * D / a2 ;            
		a1 = alfa / a ;
		MorzCoeff = 2 * a1 * D  ; 
		Ka = K * r ;                       
		K2a2 = Ka*Ka ;        
		
		for (var i = 0; i < j; i++){               
	   
				addNewC( w* 0.25 , h*0.2*(1+i), 10 );
				addNewO1(w* 0.25 + a , h*0.2*(1+i) , 10 );
				addNewO2(w* 0.25 - a, h*0.2*(1+i) , 10 );  
		  
				addNewC( w* 0.75 , h*0.2*(1+i), 10 );
				addNewO1(w* 0.75 + a , h*0.2*(1+i) , 10 );
				addNewO2(w* 0.75 - a, h*0.2*(1+i) , 10 );  
		   
		}
	}
		if(j == 20){
		j = j*0.25 ;
		r = 0.15 * a0 ;
		a = 2.5 * r;                        
	 a_1 = 2.1 * r;                       
     a_2 = 2.1 * a ;
	 aCut =  2 * r ; 
		r2 = r * r ;                      
		aCut2 = aCut * aCut ;              
		a2 = a * a ;                       
		a22 = a_2 * a_2 ;
		a11 = a_1 * a_1 ;  	
		D = a2 * Cball / 72 ;              
		LJCoeff = 12 * D / a2 ;            
		a1 = alfa / a ;
		MorzCoeff = 2 * a1 * D  ; 
		Ka = K * r ;                       
		K2a2 = Ka*Ka ;        	 
		
		for (var i = 0; i < j; i++){               
       
	   
				addNewC( w* 0.2 , h*0.16*(1+i), 10 );
				addNewO1(w* 0.2 + a , h*0.16*(1+i) , 10 );
				addNewO2(w* 0.2 - a, h*0.16*(1+i) , 10 );  
		  
				addNewC( w* 0.4 , h*0.16*(1+i), 10 );
				addNewO1(w* 0.4 + a , h*0.16*(1+i) , 10 );
				addNewO2(w* 0.4 - a, h*0.16*(1+i) , 10 );  
				
				addNewC( w* 0.6 , h*0.16*(1+i), 10 );
				addNewO1(w* 0.6 + a , h*0.16*(1+i) , 10 );
				addNewO2(w* 0.6 - a, h*0.16*(1+i) , 10 ); 
				
				addNewC( w* 0.8 , h*0.16*(1+i), 10 );
				addNewO1(w* 0.8 + a , h*0.16*(1+i) , 10 );
				addNewO2(w* 0.8 - a, h*0.16*(1+i) , 10 ); 
		   
		}
	}
		if(j == 40){
		r = 0.1 * a0 ;
		 a = 3.1 * r;                       
		a_1 = 2.1 * r;                       
		a_2 = 2.1 * a ;
		aCut =  2 * r ;  
		
		j = j*0.2 ;
		
		r2 = r * r ;                      
		aCut2 = aCut * aCut ;              
		a2 = a * a ;                       
		a22 = a_2 * a_2 ;
		a11 = a_1 * a_1 ;  	
		D = a2 * Cball / 72 ;              
		LJCoeff = 12 * D / a2 ;            
		a1 = alfa / a ;
		MorzCoeff = 2 * a1 * D  ; 
		Ka = K * r ;                       
		K2a2 = Ka*Ka ;        
		
		for (var i = 0; i < j; i++){               
				addNewC( w* 0.16 , h*0.11*(1+i), 10 );
				addNewO1(w* 0.16 + a , h*0.11*(1+i) , 10 );
				addNewO2(w* 0.16 - a, h*0.11*(1+i) , 10 );  
		  
				addNewC( w* 0.33 , h*0.11*(1+i), 10 );
				addNewO1(w* 0.33 + a , h*0.11*(1+i) , 10 );
				addNewO2(w* 0.33 - a, h*0.11*(1+i) , 10 );  
				
				addNewC( w* 0.50 , h*0.11*(1+i), 10 );
				addNewO1(w* 0.50 + a , h*0.11*(1+i) , 10 );
				addNewO2(w* 0.50 - a, h*0.11*(1+i) , 10 ); 
				
				addNewC( w* 0.67 , h*0.11*(1+i), 10 );
				addNewO1(w* 0.67 + a , h*0.11*(1+i) , 10 );
				addNewO2(w* 0.67 - a, h*0.11*(1+i) , 10 ); 
				
				addNewC( w* 0.84 , h*0.11*(1+i), 10 );
				addNewO1(w* 0.84 + a , h*0.11*(1+i) , 10 );
				addNewO2(w* 0.84 - a, h*0.11*(1+i) , 10 ); 
		   
		}
	}
		if(j == 70){
		r = 0.08 * a0 ;
		 a = 2.3 * r;                       
		a_1 = 2.1 * r;                       
		a_2 = 2.1 * a ;
		aCut =  2 * r ;  
		
		j = j/7 ;
		
		r2 = r * r ;                      
		aCut2 = aCut * aCut ;              
		a2 = a * a ;                       
		a22 = a_2 * a_2 ;
		a11 = a_1 * a_1 ;  	
		D = a2 * Cball / 72 ;              
		LJCoeff = 12 * D / a2 ;            
		a1 = alfa / a ;
		MorzCoeff = 2 * a1 * D  ; 
		Ka = K * r ;                       
		K2a2 = Ka*Ka ;        
		
		for (var i = 0; i < j; i++){               
       
				addNewC( w* 0.125 , h*0.09*(1+i), 10 );
				addNewO1(w* 0.125 + a , h*0.09*(1+i) , 10 );
				addNewO2(w* 0.125 - a, h*0.09*(1+i) , 10 );  
		  
				addNewC( w* 0.25 , h*0.09*(1+i), 10 );
				addNewO1(w* 0.25 + a , h*0.09*(1+i) , 10 );
				addNewO2(w* 0.25 - a, h*0.09*(1+i) , 10 );  
				
				addNewC( w* 0.375 , h*0.09*(1+i), 10 );
				addNewO1(w* 0.375 + a , h*0.09*(1+i) , 10 );
				addNewO2(w* 0.375 - a, h*0.09*(1+i) , 10 ); 
				
				addNewC( w* 0.5 , h*0.09*(1+i), 10 );
				addNewO1(w* 0.5 + a , h*0.09*(1+i) , 10 );
				addNewO2(w* 0.5 - a, h*0.09*(1+i) , 10 ); 
				
				addNewC( w* 0.625 , h*0.09*(1+i), 10 );
				addNewO1(w* 0.625 + a , h*0.09*(1+i) , 10 );
				addNewO2(w* 0.625 - a, h*0.09*(1+i) , 10 ); 
				
				addNewC( w* 0.75 , h*0.09*(1+i), 10 );
				addNewO1(w* 0.75 + a , h*0.09*(1+i) , 10 );
				addNewO2(w* 0.75 - a, h*0.09*(1+i) , 10 ); 
				
				addNewC( w* 0.875 , h*0.09*(1+i), 10 );
				addNewO1(w* 0.875 + a , h*0.09*(1+i) , 10 );
				addNewO2(w* 0.875 - a, h*0.09*(1+i) , 10 ); 
				}
		}
		for (var i = 0; i < C.length; i++) 
				{
					C[i].vx = 0.4 * (1 - 2 * Math.random());
                    C[i].vy = 0.4 * (1 - 2 * Math.random());
					O1[i].vx = 0.3 * (1 - 2 * Math.random());
                    O1[i].vy = 0.3 * (1 - 2 * Math.random());
					O2[i].vx = 0.3 * (1 - 2 * Math.random());
                    O2[i].vy = 0.3 * (1 - 2 * Math.random());
				}
		timers();
    }

    document.onmouseup = function(e) {          // функция при отпускании клавиши мыши
        canvas.onmousemove = null;              // когда клавиша отпущена - функции перемещения нету
        dNd = null;                             // когда клавиша отпущена - захваченного курсором шара нету
    };

    function mouseMove(e) {                     // функция при перемещении мыши, работает только с зажатой ЛКМ
        var m = mouseCoords(e);                 // получаем расчетные координаты курсора мыши
        dNd.x = m.x + dNd.xPlus;
        dNd.y = m.y + dNd.yPlus;
        dNd.vx = 0.6 * (m.x - mx_) / dt / fps;   dNd.vy = 0.6 * (m.y - my_) / dt / fps;
        mx_ = m.x;    my_ = m.y;
    }

    function mouseCoords(e) {                   // функция возвращает расчетные координаты курсора мыши
        var m = [];
        var rect = canvas.getBoundingClientRect();
        m.x = (e.clientX - rect.left) / scale;
        m.y = (e.clientY - rect.top) / scale;
        return m;
    }

    

	var addNewBall =  function(x, y) {
        // проверка - не пересекается ли новый шар со стенами или уже существующими шарами
        if (x - r < 0 || x + r > w || y - r < 0 || y + r > h) return null;
        for (var i = 0; i < balls.length; i++) {
            var rx = balls[i].x - x;
            var ry = balls[i].y - y;
            var rLen2 = rx * rx + ry * ry;
            if (rLen2 < 4 * r2) return null;
        }
        var b = [];

        b.x = x;                b.y = y;        // расчетные координаты шара
        b.fx = 0;               b.fy = 0;      // сила, действующая на шар
        b.vx = 0;               b.vy = 0;   		// скорость 

        balls[balls.length] = b;                // добавить элемент в конец массива
        
		return b; 
    };
	
    var addNewC =  function(x, y) {
        // проверка - не пересекается ли новый шар со стенами или уже существующими шарами
        if (x - r < 0 || x + r > w || y - r < 0 || y + r > h) return null;
        for (var i = 0; i < balls.length; i++) {
            var rx = balls[i].x - x;
            var ry = balls[i].y - y;
            var rLen2 = rx * rx + ry * ry;
            if (rLen2 < 4 * r2) return null;
        }
        var b = [];

        b.x = x;                b.y = y;        // расчетные координаты шара
        b.fx = 0;               b.fy = 0;      // сила, действующая на шар
        b.vx = 0;               b.vy = 0;   		// скорость

        C[C.length] = b;                // добавить элемент в конец массива
        
		return b; 
    };

    var addNewO1 =  function(x, y) {
        // проверка - не пересекается ли новый шар со стенами или уже существующими шарами
        if (x - r < 0 || x + r > w || y - r < 0 || y + r > h) return null;
        for (var i = 0; i < balls.length; i++) {
            var rx = balls[i].x - x;
            var ry = balls[i].y - y;
            var rLen2 = rx * rx + ry * ry;
            if (rLen2 < 4 * r2) return null;
        }
        var b = [];

        b.x = x;                b.y = y;        // расчетные координаты шара
        b.fx = 0;               b.fy = 0;      // сила, действующая на шар
        b.vx = 0;               b.vy = 0;   		// скорость
		
        O1[O1.length] = b;                // добавить элемент в конец массива
        
		return b; 
    };
	
	var addNewO2 =  function(x, y) {
        // проверка - не пересекается ли новый шар со стенами или уже существующими шарами
        if (x - r < 0 || x + r > w || y - r < 0 || y + r > h) return null;
        for (var i = 0; i < balls.length; i++) {
            var rx = balls[i].x - x;
            var ry = balls[i].y - y;
            var rLen2 = rx * rx + ry * ry;
            if (rLen2 < 4 * r2) return null;
        }
        var b = [];

        b.x = x;                b.y = y;        // расчетные координаты шара
        b.fx = 0;               b.fy = 0;      // сила, действующая на шар
        b.vx = 0;               b.vy = 0;   		// скорость

        O2[O2.length] = b;                // добавить элемент в конец массива
        
		return b; 
    };
	
	// Основной цикл программы

    function control() {
        physics();
        draw();
    }
	function timers(){
		clearInterval(id1);
		clearInterval(id2);
		clearInterval(id3);
		clearInterval(id4);
		clearInterval(id5);
		k = 1;
		Ekc = 0 ;
		Ekck = 0;
		EC= 0;
		EO1 = 0 ;
		EO2 = 0;
		ECk = 0;
		EO1k =0 ; 
		EO2k = 0 ;
		Ek = 0;
		Ek1 = 0 ;
		
		var id1 = setInterval(function(){document.getElementById('Ekck').innerHTML = Ekck;}, 1000 / 2);
		var id2 =setInterval(function(){document.getElementById('Ek').innerHTML = Ek;}, 1000 / 2);
		var id3 =setInterval(function(){document.getElementById('Ekc').innerHTML = ECk;}, 1000 / 2);
		var id4 =setInterval(function(){document.getElementById('Eko1').innerHTML = EO1k;}, 1000 / 2);
		var id5 =setInterval(function(){document.getElementById('Eko2').innerHTML = EO2k;}, 1000 / 2);
		
	}
    // Расчетная часть программы

	function powers(b ,b2 , hkk , jk) {
		var rx = b.x - b2.x;   var ry = b.y - b2.y;         // вектор смотрит на первый шар (b)
        var r2 = rx * rx + ry * ry;                         // квадрат расстояния между шарами
                  	
					if ( r2 > aCut2 ) return null;	
					
					var rLen = (Math.sqrt(r2));
       
					if (r2 < K2a2) {
                        if (rLen > 0.00001) {                           // проверка, чтобы избежать деления на 0
                            rx = rx / rLen * Ka;
                            ry = ry / rLen * Ka;
                        }
                        r2 = K2a2;
                        rLen = Ka;                                      // корень K2a2
                    }

                    // сила взаимодействия
                   
					var u = Math.exp( -a1 * ( rLen - hkk  )) ;
					var F = jk  * MorzCoeff * u* ( u - 1 ) /rLen ;
					
                    var Fx = F * rx;        var Fy = F * ry;
                    b.fx += Fx;             b.fy += Fy;
                    b2.fx -= Fx;            b2.fy -= Fy;	
	}
	
	
    function physics() {                        // то, что происходит каждый шаг времени
        for (var s = 1; s <= spf; s++) {

            // пересчет сил идет отдельным массивом, т.к. далее будут добавляться силы взаимодействия между шарами
           
		   for (var i0 = 0; i0 < balls.length; i0++) {
                balls[i0].fx = - B * balls[i0].vx;
                balls[i0].fy =  - B * balls[i0].vy;
            }
			for (var i0 = 0; i0 < C.length; i0++) {
                C[i0].fx = - B * C[i0].vx;
                C[i0].fy =  - B * C[i0].vy;
            }
			for (var i0 = 0; i0 < O1.length; i0++) {
                O1[i0].fx = - B * O1[i0].vx;
                O1[i0].fy =  - B * O1[i0].vy;
            }
			for (var i0 = 0; i0 < O2.length; i0++) {
                O2[i0].fx = - B * O2[i0].vx;
                O2[i0].fy =  - B * O2[i0].vy;
            }

			
            for (var i = 0; i < balls.length; i++) { // dlya Balls
                
				var b = balls[i];
                
				for (var j = i + 1; j < balls.length; j++) {
                    var b2 = balls[j];
					powers(b,b2,a_1,0.3) ;
				}
				 
				 for (var j = 0; j < C.length; j++) {
                    var b2 = C[j];
					powers(b,b2,a_1,0.3) ;
				
				}
				 for (var j = 0; j < O1.length; j++) {
                    var b2 = O1[j];
					powers(b,b2,a_1,0.3) ; 
				}
				
				 for (var j = 0; j < O2.length; j++) {
                    var b2 = O2[j];
					powers( b, b2, a_1 ,0.3) ;
				}
				
				if (b == dNd) continue;  // если шар схвачен курсором - его вз. со стенами и перемещение не считаем

                if (b.y + r > h) { b.fy += -Cwall * (b.y + r - h) - Bwall * b.vy; }
                if (b.y - r < 0) { b.fy += -Cwall * (b.y - r) - Bwall * b.vy;}
                if (b.x + r > w) { b.fx += -Cwall * (b.x + r - w) - Bwall * b.vx; }
                if (b.x - r < 0) { b.fx += -Cwall * (b.x - r) - Bwall * b.vx; }

                b.vx += b.fx / m * dt;        b.vy += b.fy / m * dt;
                b.x += b.vx * dt;             b.y += b.vy * dt;
				Ek1 += m1* 0.5 *(b.vx*b.vx + b.vy*b.vy);
            }
			
			for (var i = 0; i < C.length; i++) {
				for (var j = i ; j < C.length; j++) {
					
					if ( i == j){
						var r1x = O1[i].x - C[i].x ; var r1y = O1[i].y - C[i].y ;
						var r2x = O2[i].x - C[i].x ; var r2y = O2[i].y - C[i].y ; 
					
						var r1 = Math.sqrt(r1x*r1x +r1y* r1y);
						var r2 = Math.sqrt(r2x*r2x +r2y* r2y);
						var xx =(r1x*r2x + r1y*r2y)/(r1*r2) ;
					
						var F1 = -0.2*Cball*(r1-a)/r1;
						var F2 = -0.2*Cball*(r2-a)/r2;
						var f = 0 ;
						if(xx*xx > 0.99999999999){
							f = -0.3*Cball;	
							}
						
						else{
						f = -0.3*Cball*(Pi - Math.acos(xx))/(Math.sqrt(1-xx*xx));
						}
						ff = f ;
						
						var f1x = f*(r2x/(r1*r2) - xx*r1x/(r1*r1)); var f1y = f*(r2y/(r1*r2) -xx*r1y/(r1*r1));
						var f2x = f*(r1x/(r1*r2) - xx*r2x/(r2*r2)); var f2y = f*(r1y/(r1*r2) -xx*r2y/(r2*r2));
					
						var F1x = F1 * r1x   ;        var F1y = F1 * r1y ;
						var F2x = F2 * r2x   ;        var F2y = F2 * r2y ;
						
						O1[i].fx += (F1x + f1x) ;  O1[i].fy += (F1y + f1y) ;
						O2[i].fx += (F2x + f2x) ;  O2[i].fy += (F2y + f2y) ; 
						C[i].fx -= (F1x + F2x + f1x + f2x) ;  C[i].fy -= (F1y + F2y + f1y + f2y) ; 
						
						powers(C[i],C[j],a_1, 0.3) ; powers(C[i],O1[j],a_1, 0.3) ; powers(C[i],O2[j],a_1, 0.3) ;
					powers(O1[i],C[j],a_1, 0.3) ; powers(O1[i],O1[j],a_1, 0.3) ; powers(O1[i],O2[j],a_1, 0.3) ;
					powers(O2[i],C[j],a_1, 0.3) ; powers(O2[i],O1[j],a_1, 0.3) ; powers(O2[i],O2[j],a_1, 0.3) ;
						
						}
					
					powers(C[i],C[j],a_1, 0.3) ; powers(C[i],O1[j],a_1, 0.3) ; powers(C[i],O2[j],a_1, 0.3) ;
					powers(O1[i],C[j],a_1, 0.3) ; powers(O1[i],O1[j],a_1, 0.3) ; powers(O1[i],O2[j],a_1, 0.3) ;
					powers(O2[i],C[j],a_1, 0.3) ; powers(O2[i],O1[j],a_1, 0.3) ; powers(O2[i],O2[j],a_1, 0.3) ;
					
				}
			}
			
			
			
			for (var i = 0; i < C.length; i++) { // dlya C
                 
				var b = C[i];
                
				if (b == dNd) continue;  // если шар схвачен курсором - его вз. со стенами и перемещение не считаем

                if (b.y + r > h) { b.fy += -Cwall * (b.y + r - h) - Bwall * b.vy; }
                if (b.y - r < 0) { b.fy += -Cwall * (b.y - r) - Bwall * b.vy;}
                if (b.x + r > w) { b.fx += -Cwall * (b.x + r - w) - Bwall * b.vx; }
                if (b.x - r < 0) { b.fx += -Cwall * (b.x - r) - Bwall * b.vx; }

             
                b.vx += b.fx / m * dt;        b.vy += b.fy / m * dt;
                b.x += b.vx * dt;             b.y += b.vy * dt;
				Ek1 += m* 0.5 *(b.vx*b.vx + b.vy*b.vy) ;
				EC += m* 0.5 *(b.vx*b.vx + b.vy*b.vy)/j ;
            }
			
			for (var i = 0; i < C.length; i++) { // dlya O1
              
				var b = O1[i];
     
				if (b == dNd) continue;  // если шар схвачен курсором - его вз. со стенами и перемещение не считаем

                if (b.y + r > h) { b.fy += -Cwall * (b.y + r - h) - Bwall * b.vy; }
                if (b.y - r < 0) { b.fy += -Cwall * (b.y - r) - Bwall * b.vy;}
                if (b.x + r > w) { b.fx += -Cwall * (b.x + r - w) - Bwall * b.vx; }
                if (b.x - r < 0) { b.fx += -Cwall * (b.x - r) - Bwall * b.vx; }

                b.vx += b.fx / m1 * dt;        b.vy += b.fy / m1 * dt;
                b.x += b.vx * dt;             b.y += b.vy * dt;
				Ek1 += m1 * 0.5 *(b.vx*b.vx + b.vy * b.vy) ;
				EO1 += m1* 0.5 *(b.vx*b.vx + b.vy*b.vy)/j ;
			
                
            }
			
			for (var i = 0; i < C.length; i++) { // dlya O2
                
				var b = O2[i];
                
				 
				if (b == dNd) continue;  // если шар схвачен курсором - его вз. со стенами и перемещение не считаем

                if (b.y + r > h) { b.fy += -Cwall * (b.y + r - h) - Bwall * b.vy; }
                if (b.y - r < 0) { b.fy += -Cwall * (b.y - r) - Bwall * b.vy;}
                if (b.x + r > w) { b.fx += -Cwall * (b.x + r - w) - Bwall * b.vx; }
                if (b.x - r < 0) { b.fx += -Cwall * (b.x - r) - Bwall * b.vx; }

                b.vx += b.fx / m1 * dt;        b.vy += b.fy / m1 * dt;
                b.x += b.vx * dt;             b.y += b.vy * dt;
				Ek1 += m1* 0.5 *(b.vx*b.vx + b.vy*b.vy) ;
				EO2 += m1* 0.5 *(b.vx*b.vx + b.vy*b.vy)/j ;
				
				
            
			}
			 
        }
		Ekc +=Ek1 ;
		
		if(Ek1 > 1){
			k++;
			}
		Ekck = Ekc/k;
		ECk = EC/k ;
		EO1k = EO1/k ; 
		EO2k = EO2/k ;
		Ek = Ek1;
		Ek1 = 0 ;
    }
	
   
   

   // Рисование

    var rScale13 = r * scale * 1.3;         // ___в целях оптимизации___
    var rScaleShift = r * scale / 5;        // ___в целях оптимизации___
    function draw() {
        context.clearRect(0, 0, w * scale, h * scale);      // очистить экран
        for (var i = 0; i < balls.length; i++){
            var xS = balls[i].x * scale;           var yS = balls[i].y * scale;
            
            context.fillStyle = "#FFD700";

            context.beginPath();
            context.arc(xS, yS, r * scale, 0, 2 * Math.PI, false);
            context.closePath();
            context.fill();
        }
    
		for (var i = 0; i < C.length; i++){
            var xS = C[i].x * scale;           var yS = C[i].y * scale;
            
            context.fillStyle = "#f08080";

            context.beginPath();
            context.arc(xS, yS, r * scale, 0, 2 * Math.PI, false);
            context.closePath();
            context.fill();
        }
		for (var i = 0; i < O1.length; i++){
            var xS = O1[i].x * scale;           var yS = O1[i].y * scale;
            
            context.fillStyle = "#3070d0";

            context.beginPath();
            context.arc(xS, yS, r * scale, 0, 2 * Math.PI, false);
            context.closePath();
            context.fill();
        }
		for (var i = 0; i < O2.length; i++){
            var xS = O2[i].x * scale;           var yS = O2[i].y * scale;
            
            context.fillStyle = "#3070d0";

            context.beginPath();
            context.arc(xS, yS, r * scale, 0, 2 * Math.PI, false);
            context.closePath();
            context.fill();
        }
	}

    // Запуск системы
	
	this.newSystem();
	var id1 = setInterval(function(){document.getElementById('Ekck').innerHTML = Ekck;}, 1000 / 2);
	var id2 =setInterval(function(){document.getElementById('Ek').innerHTML = Ek;}, 1000 / 2);
	var id3 =setInterval(function(){document.getElementById('Ekc').innerHTML = ECk;}, 1000 / 2);
	var id4 =setInterval(function(){document.getElementById('Eko1').innerHTML = EO1k;}, 1000 / 2);
	var id5 =setInterval(function(){document.getElementById('Eko2').innerHTML = EO2k;}, 1000 / 2);
	
    setInterval(control, 1000 / fps);
}
</syntaxhighlight>
Файл '''"Lab.html"'''
<syntaxhighlight lang="html5" line start="1" enclose="div">
<!DOCTYPE html>
<html>
<head>
    <title>Balls</title>
    <script src="Balls_v4_release.js"></script>
</head>
<body>
    <canvas id="canvasBalls" width="640" height="480" style="border:1px solid #000000;"></canvas>
	<br>
	 <div>Количество частиц:
        <input type="button" style="width: 30px" name="" onclick="setNy(8); newSystem();return false;" value="8"/>
        <input type="button" style="width: 30px" name="" onclick="setNy(20); newSystem();return false;" value="20"/>
        <input type="button" style="width: 30px" name="" onclick="setNy(40); newSystem();return false;" value="40"/>
        <input type="button" style="width: 30px" name="" onclick="setNy(70); newSystem();return false;" value="70"/>
    </div>
	 <div>Обновить таймер:
        <input type="button" style="width: 60px" name="" onclick="ResetTime();return false;" value="Reset"/>
    </div>
	<div>Кин энергия в данный момент времени: <span id="Ek"></span></div>
	<div>Средняя кинетическая энергия: <span id="Ekck"></span></div>
	<div>Средняя кинетическая энергия на один атом углерода: <span id="Ekc"></span></div>
	<div>Средняя кинетическая энергия на один атом кислорода(1): <span id="Eko1"></span></div>
	<div>Средняя кинетическая энергия на один атом кислорода(2): <span id="Eko2"></span></div>
</body>
</html>
</syntaxhighlight>
</div>
</div>

==Результаты==
Получена программа, которая моделирует молекулы углекислого газа в объеме и рассчитывает среднюю кинетическую энергию отдельных частей системы и всей системы в целом.
Значение кинетической энергии системы совершает колебания вокруг среднего значения энергии, которое оказалось постоянным, что говорит о выполнении закона сохранения энергии. Средние энергии, приходящиеся на каждый отдельный атом, спустя какое-то время оказываются равными. Это говорит о том, что энергии, проходящиеся на каждую степень свободы, равны. Для идеального газа эта энергия равна:
::<math>U= \frac{kT}{2} </math>
где
* <math>k</math> — постоянная Больцмана,
* <math>T</math> — температура среды,  
А нашу систему можно приблизительно считать идеальным газом в силу малости значения потенциальной энергии и абсолютно упругих соударений. Значит наши результаты соответствуют, теории.