栏目分类:
子分类:
返回
名师互学网用户登录
快速导航关闭
当前搜索
当前分类
子分类
实用工具
热门搜索
名师互学网 > IT > 软件开发 > 后端开发 > C/C++/C#

[C++调试笔记]推动粒子move.cpp

C/C++/C# 更新时间: 发布时间: IT归档 最新发布 模块sitemap 名妆网 法律咨询 聚返吧 英语巴士网 伯小乐 网商动力

[C++调试笔记]推动粒子move.cpp

#include 
#include 
#include "define.h"
void move(particle* ptr, int* numb, double* fiel, double q, double m, double* flux, double* ener_flux)
{
	int i, j, ip;
	double tt[3], u[3], tt1, u1[3], u2[3], u3[3], tan_b;	//存放Boris算法中的各速度分量信息
	double z, ez, az;
	double l1, l2;
	double vx, vy, vz;
	double weig, ener;
	double cos_thet;							//入射方向与靶板法线方向夹角cos值

	j = 0;										//统计还在计算区域内的粒子数
	for (i = 0; i < 2; i++)
	{
		flux[i] = 0.0;
		ener_flux[i] = 0.0;
	}

	for (i = 0; i < *numb; i++)
	{
		z = ptr[i].z;
		vx = ptr[i].vx;
		vy = ptr[i].vy;
		vz = ptr[i].vz;
		weig = ptr[i].weig;
		
		ip = (int)(z / dz_plas);
		l2 = z / dz_plas - ip;
		l1 = 1.0 - l2;
		ez = fiel[ip] * l1 + fiel[ip + 1] * l2;
		az = q * ez / m;
		
		vz = vz + az * dt;

		z = z + vz * dt;

		
		if ((z > 0.0) && (z < lz_plas))
		{
			ptr[j].z = z;
			ptr[j].vx = vx;
			ptr[j].vy = vy;
			ptr[j].vz = vz;
			ener = 0.5 * m * (vx * vx + vy * vy + vz * vz) / qi;
			ptr[j].ener = ener;
			ptr[j].weig = weig;
			j++;
		}
		else if (z <= 0.0)
		{
			ener = 0.5 * m * (vx * vx + vy * vy + vz * vz) / qi;												//单位eV,,需要进行转换
			cos_thet = fabs(vz) / sqrt(vx * vx + vy * vy + vz * vz);
			flux[0] = (flux[0] + weig) * cos_thet;
			ener_flux[0] = (ener_flux[0] + weig * ener) * cos_thet;

		}
		else if (z >= lz_plas)
		{
			ener = 0.5 * m * (vx * vx + vy * vy + vz * vz) / qi;												//单位eV,,需要进行转换
			cos_thet = fabs(vz) / sqrt(vx * vx + vy * vy + vz * vz);
			flux[1] = (flux[1] + weig) * cos_thet;
			ener_flux[1] = (ener_flux[1] + weig * ener) * cos_thet;

		}
	}

	*numb = j;
	for (i = 0; i < 2; i++)
	{
		flux[i] = flux[i] / (area * dt);
		ener_flux[i] = ener_flux[i] * qi / (area * dt);
	}
}

转载请注明:文章转载自 www.mshxw.com
本文地址:https://www.mshxw.com/it/347028.html
我们一直用心在做
关于我们 文章归档 网站地图 联系我们

版权所有 (c)2021-2022 MSHXW.COM

ICP备案号:晋ICP备2021003244-6号