matlab如何做粒子模拟,求助,如何用matlab做蒙特卡罗模拟!!??
我这儿正好有份程序,希望有所帮助。
代码:
一份蒙特卡洛程序
count=input('input the count:');%输入模拟粒子数
sigmaedata=[28370,13845,6908,2555,1223,2602,1925,905.3,479.7,164.5,74.24,23.86,66.60,36.62,22.29,9.978,5.298,1.668,0.7378,0.2361,0.1099,0.06211,0.03939,0.02030];
sigmacdata=[0.0220,0.0393,0.0568,0.0904,0.1209,0.1479,0.1722,0.2136,0.2480,0.3092,0.3486,0.3932,0.4153,0.4268,0.4319,0.4291,0.4215,0.3969,0.3691,0.3269,0.2944,0.2709,0.2512,0.2209];
Edata=[1,1.5,2,3,4,5,6,8,10,15,20,30,40,50,60,80,100,150,200,300,400,500,600,800];%截面数据
channel=zeros(1,ceil(662/5)+10);%多道数组
nget=0;%探测到的总计数
ntotal=0;%进入探测器的总计数
for ii=1:count%count个粒子循环
collidetime=0;%当前粒子碰撞次数
%粒子状态初始化
E0=622;
E=E0;
z=-2;
r=0;
theta=2*pi*rand(1);% 源抽样,z,r,theta坐标
miu=2*rand(1)-1;
fai=2*pi*rand(1);%方向角抽样
if miu
continue;
else
z=0;
r=2*sqrt(1-miu^2)/miu;
theta=fai;
end
while E>1%一个粒子在闪烁体中的输运过程
sigmae=interp1(Edata,sigmaedata,E,'linear');
sigmac=interp1(Edata,sigmacdata,E,'linear');
sigmat=sigmae+sigmac;%线性插值得到截面数据
L=-log(rand(1))/sigmat;%下次碰撞的距离
%计算下次碰撞位置坐标
rnew=sqrt(r^2+L^2*(1-miu^2)+2*r*L*sqrt(1-miu^2)*cos(fai-theta));
z=z+L*miu;
cdth=(rnew^2+r^2-L^2*(1-miu^2))/2/r/rnew;
sdth=L*sqrt(1-miu^2)*sin(fai-theta)/rnew;
dtheta=asin(sdth);
if cdth<0
dtheta=pi-dtheta;
end
theta=theta+dtheta;
r=rnew;
if(r>2)|(z>=4)|(z<0)%判断是否在闪烁体内
break;
else
collidetime=collidetime+1;
end
if rand(1)
E=0;
else%康普顿散射
alpha=E/511;
flag=0;
while flag==0
if rand(1)<=27/(4*alpha+29)
x=(1+2*alpha)/(1+2*alpha*rand(1));
if rand(1)<=0.5*((alpha+1-x/alpha)^2+1)
flag=1;
end
else
x=1+2*alpha*rand(1);
if rand(1)<=27/4*((x-1)^2)/x^3
flag=1;
end
end
end
E=E/x;
alphat=alpha/x;
miuL=1-1/alphat+1/alpha;%散射后的方向
a=miuL;
b=sqrt(1-a^2);
randangle=2*pi*rand(1);
miunew=a*miu+b*sqrt(1-miu^2)*cos(randangle);
sdf=b*sin(randangle)/sqrt(1-miunew^2);
cdf=(a-miu*miunew)/sqrt(1-miu^2)/sqrt(1-miunew^2);
sfn=sdf*cos(fai)+cdf*sin(fai);
cfn=cdf*cos(fai)-sdf*sin(fai);
fainew=asin(sfn);
if(cfn<0)
fainew=pi-fainew;
end
fai=fainew;
miu=miunew;
end%conputon/electtron
%记录结果
end%while
ntotal=ntotal+1;%进入探测器的总计数
if collidetime>0
nget=nget+1;%探测到的计数
Edown=E0-E;%沉积能量
FWHM=0.01+0.05*sqrt(Edown+0.4*Edown^2);
delta=0.4247*FWHM;
Eget=Edown+delta*randn(1);%记录能量
num=ceil(Eget/5); %寻道址
if num~=0
channel(num)=channel(num)+1;
end
end
end%for count
plot(channel);
fprintf(1,'探测效率为%1.5f',nget/ntotal);
matlab如何做粒子模拟,求助,如何用matlab做蒙特卡罗模拟!!??相关推荐
- matlab引用数据,excel引用数据-如何用matlab处理excel文件中的数据?
如何利用matlab根据excel表格里面的数据画图 将待导入的矩阵结数据Excel中,录入时注意行列原矩阵一一对应 录入完以后数据,为了后续步骤使用方便,命名时我们最好把它命名为我们接下来在MATL ...
- matlab鼠标三维坐标点,请问如何用matlab画三维点,已知x,y,z的坐标,在三维坐标系上显示...
点击查看请问如何用matlab画三维点,已知x,y,z的坐标,在三维坐标系上显示具体信息 答:例如 : X=1,Y=2,Z=3; 代码就是: plot3(1,2,3,'*') grid on%加网格 ...
- matlab求分段函数的值.,如何用MATLAB求分段函数的最小值和最大值?
7.1.1 分段线性插值 所谓分段线性插值就是通过插值点用折线段连接起来逼近原曲线,这也是计算机绘制图形的基本原理.实现分段线性插值不需编制函数程序,MATLAB自身提供了内部函数interp1其主要 ...
- PR做片头:如何用PR做片头,如何用模板来制作片头?
如何用PR做片头?本站提供了海量的PR模板,学会利用模板来制作,简单又高效,那么给大家讲讲PR做片头过程中的一些小技巧. 先在本网站上下载需要的PR模板,下载到电脑之后解压 PR模板:片头片尾 打开工 ...
- 如何用R进行蒙特卡罗模拟(Monte Carlo Simulation with R)
本文所讲的蒙特卡罗模拟是建立在正态分布的基础上.假设我们给定一只股票的初始价格P0, 并且从历史日度数据估计出该股票的每日期望回报率为 mean.logret, 标准差为sd.logret. ( 在估 ...
- matlab模块里有s,求助!!S-Function做通用模块
三.简单单摆S-函数仿真实验 使用S-函数来对一个单摆系统进行仿真,主要演示以下三个方面: (1)利用S-函数对单摆系统进行建模 (2)利用simulink进行仿真,研究单摆的位移 (3)利用S-函数 ...
- 代入消元法 matlab,求助 如何用matlab计算期权价格
期权定价理论是现代金融学中最为重要的理论之一,也是衍生金融工具定价中最复杂的.本文给出了欧式期权定价过程的一个简单推导,并利用Matlab对定价公式给出了数值算例及比较静态分析,以使读者能更直观地理解 ...
- matlab 怎么与运算,求助如何用MATLAB计算VAR和Expected shortfall
不知道这个能不能满足 function calcvar %calcvar terminal interface for VaR calculation. % % You c ...
- matlab找距离最近的元素,如何用MATLAB找到给定坐标的最近点?
这总是很有趣:) 首先:Mohsen Nosratinia的回答是好的,只要 您不需要知道实际距离 你可以绝对肯定地保证你永远不会去极地附近 并且永远不会接近±180°子午线 对于给定的纬度,-180 ...
- matlab 图案 柱状图_值得收藏 | 如何用matlab做出酷炫的图像
1.基础知识 1.1 二维图形绘制 plot函数是Matlab绘制二维图形的常用函数,该函数将数组中的数据点绘连起来构成一条连续的曲线. plot(x,y,'PropertyName',Propert ...
最新文章
- 注入游戏没有焦点_不戴眼镜看3D电影、玩3D游戏,这项技术能焕发端游市场第二春吗?...
- Linux网络模拟,模拟网络访问解析
- nyoj655光棍的yy
- web.py——运行错误【AttributeError: ‘StaticApp‘ object has no attribute ‘directory‘】
- 模块开发卷宗(GB8567——88)
- halcon write_ocr_class_svm 将OCR分类器写入文件
- Linux卸载MariaDB
- git将本地仓库推送到远程仓库
- 魅族MX4关闭系统升级Flyme6提示
- 为什么div设置其border无效?
- SQLServer 2008 r2 下载地址(百度云)及安装图解
- 使用wbadmin备份整个网络上的完整Vista PC
- 如何解决merge conflict的方法
- Linux 学习包括但不限于linux使用问题笔记
- Android使用App Architecture打造最佳体验和高质量应用《一》
- python-matplotlib 绘制函数曲线
- 拓新药业301089
- 1024程序员节的由来
- 想在美国创业却没有H1B?这些大学可以帮忙搞定身份!
- 如何从CCleaner清理中排除项目
热门文章
- Java复习题1(1.写出抽象类和接口的区别。2.This和Super关键字的区别有哪些?3.常见的类集子类有哪些,各有什么特点?4.什么是多态,多态实现的前提是什么?)
- Java唐诗学习系统
- GifCam2.0使用
- Cognos 11.0快速开发指南 Ⅰ
- 陈建文资料介绍:与众不同的交易机遇
- [Android] ListView实现隔行变色(一)
- 恶意软件与反病毒技术
- Eclipse Shell for Plugin
- python爬空气污染实时数据_一键爬取空气质量相关指数
- Spark 场景题详解