作者mr24 (=大头=)
看板Fortran
标题[请益] 程式问题请教
时间Wed Apr 25 18:23:30 2012
先PO一下我原做这程式的动机
我是有一个档案他会记录一个碳球由60颗碳所组成的空间座标档
他会随着我跑模拟後而计录每一步阶的原子位置变化
而我的目的是去计算这颗碳球的表面积，他是由20个六角形跟12个五角形所组成，
我已经将角形是由哪几颗原子所组成的纪录好，只差带入公式去计算，而我在这是采用
海龙公式。
这边开始PO我的主要程式部分，我只先PO其中一小块面积的计算部分顺便说明参数定义
iTimeStep 表我总共跑了多少步阶
iBins 表我有多少颗原子
n 是为了先给我的原子档案上编号
x,y,z我是替每一颗原子所上的编号
open(33,file='input.dat',status='unknown')
open(66,file='out.dat',status='unknown')
open(67,file='area.dat',status='unknown')
n=0
iTimeStep=0
321 format(1x, i3,1x, a1, 3(1x,e15.8))
322 format(1x,i6,1x,e15.8)
do while (1 .eq. 1)
read(33,*,end=25) iBins
read(33,*)
iTimeStep=iTimeStep+1
write (66,*)iTimeStep
write (66,*)iBins
do i = 1, iBins
read(33,*) c,x(i),y(i),z(i)
n=n+1
write(66,321)n,c,x(i),y(i),z(i)
!!aa1
if(i.eq.1)then!!!!!!!!!!!!!!!!!!!!!!!!!!!A
x1=x(i)
y1=y(i)
z1=z(i)
endif
if(i.eq.2)then
x2=x(i)
y2=y(i)
z2=z(i)
endif
if(i.eq.3)then
x3=x(i)
y3=y(i)
z3=z(i)
endif
if(i.eq.4)then
x4=x(i)
y4=y(i)
z4=z(i)
endif
if(i.eq.5)then
x5=x(i)
y5=y(i)
z5=z(i)
endif
里面的 if 是希望能够读原子编号，i=1时读编号一的原子座标，依此类推...
当读完时候开始做里面积算
cx=(x1+x2+x3+x4+x5)/5
cy=(y1+y2+y3+y4+y5)/5
cz=(z1+z2+z3+z4+z5)/5
!calculate the dis from center to each point
lc1=(((cx-x1)**2)+((cy-y1)**2)+((cz-z1)**2))**0.5
lc2=(((cx-x2)**2)+((cy-y2)**2)+((cz-z2)**2))**0.5
lc3=(((cx-x3)**2)+((cy-y3)**2)+((cz-z3)**2))**0.5
lc4=(((cx-x4)**2)+((cy-y4)**2)+((cz-z4)**2))**0.5
lc5=(((cx-x5)**2)+((cy-y5)**2)+((cz-z5)**2))**0.5
!calculate the dis from point to point
l12=(((x1-x2)**2)+((y1-y2)**2)+((z1-z2)**2))**0.5
l23=(((x2-x3)**2)+((y2-y3)**2)+((z2-z3)**2))**0.5
l34=(((x3-x4)**2)+((y3-y4)**2)+((z3-z4)**2))**0.5
l45=(((x4-x5)**2)+((y4-y5)**2)+((z4-z5)**2))**0.5
l51=(((x5-x1)**2)+((y5-y1)**2)+((z5-z1)**2))**0.5
s1=(lc1+lc2+l12)/2
s2=(lc2+lc3+l23)/2
s3=(lc3+lc4+l34)/2
s4=(lc4+lc5+l45)/2
s5=(lc5+lc1+l51)/2
!!!!!!!!!!!计算每一三角形面积
a1=(s1*(s1-lc1)*(s1-lc2)*(s1-l12))**0.5
a2=(s2*(s2-lc2)*(s2-lc3)*(s2-l23))**0.5
a3=(s3*(s3-lc3)*(s3-lc4)*(s3-l34))**0.5
a4=(s4*(s4-lc4)*(s4-lc5)*(s4-l45))**0.5
a5=(s5*(s5-lc5)*(s5-lc1)*(s5-l51))**0.5
!!!!!五边形面积
aa1=a1+a2+a3+a4+a5
enddo
t5a=aa1
totala=t5a
write(67,322)iTimeStep, totala
n=0
enddo
25 continue
end
===================================================================
我每一个原子都给他不同变数，但计算部分例如cx,cy,cz,lc1,l23,s2,a1
等共同部分则没有个别设立
目前单独计算一块面积时跟我用excel(一开始不会写只好用excel去计算)算时
几乎一样，小数三位後才不同，但我开始计算多块加总後却差别越来越大
有怀疑是否之後的计算面积是否真的有代入新座标，还是仍用第一组座标
还有不知我共同计算部分的参数是否也该个别独立宣告，也怀疑是否是这边
共同计算部分的问题。
内容有点多，希望板上的高手能帮忙指点一下，最近很烦恼 >"< 万分感谢!
--
※ 发信站: 批踢踢实业坊(ptt.cc)
◆ From: 140.116.75.123
1F：推 terryys:算面积的loop应该是先读取五个坐标以後再算面积？现在 04/26 11:34
2F：→ terryys:看起来是读一个坐标算一次 04/26 11:34
3F：→ terryys:还有totala那里应该是总面积吗？ 04/26 11:35
4F：→ mr24:totala是算面积没错，我也有想过他是不是真的有读完 04/26 13:15
5F：→ mr24:所以我计算面积部分是应该放在Do loop的外面吗?? 04/26 13:16
6F：→ mr24:还是应该再多加一个Do 04/26 13:16
7F：推 terryys:总面积的话你现在没有累加啊，只是一个五角形的面积，不过 04/26 13:21
8F：→ terryys:这不是大问题 04/26 13:22
9F：→ terryys:我觉得应该先读取iBins里的那几个坐标，然後再用这些坐标 04/26 13:22
10F：→ terryys:算面积 还有那堆if应该不需要，直接用x(1) x(2)之类就可以 04/26 13:24
11F：→ mr24:好 我去试试看，真的非常谢谢你的解答!! 04/26 14:04