作者geeoh (機機北北北北)
看板CodeJob
標題Fw: [請益] 拜託高手請幫我把fortran改為C++程式碼
時間Tue Nov 29 19:33:18 2011
※ [本文轉錄自 Mechanical 看板 #1ErB-yho ]
作者: geeoh (機機北北北北) 看板: Mechanical
標題: [請益] 拜託高手請幫我把fortran改為C++程式碼
時間: Tue Nov 29 19:19:22 2011
program nozzle
implicit real*8(a-h,o-z)
real k,kco,kh2o,kh2,kco2,num,11,keq,jh2
dimension x(8), ppt(8),ff(8),f(8),cp(6)
open(11,file='xco-k-1000.txt',status='unknown')
pi=3.1415926
dt=2*1.4625e-3
ds=3*1.46e-3
at=0.25*pi*dt**2
as=0.25*pi*(ds**2-dt**2)
ll=0.130175
rr=8.314
hinf=0.
pref=1.01e5
ptt=pref
c pmh2=4.4e-5
c pmh2=5.4e-5
c pmh2=3.4e-5
c pmh2=0.
pmh2=6.4e-5
pmh2=1000e-5
pts=1*pref
c pts=3*pref
c pts=5*pref
2001 vvlt=3.33e-6/4
vvls=1e-6
z=0.
hts=40.
c
c 1=co ,2=h2o, 3=h2, 4=co2, 5=he in permeate side, 6=h2 in permearate side,
c 7=temperature in tube side, 8=temperature in permeate side
c
c sc1.5
x(1)=0.295
x(2)=0.443
c
c sc=4
c x(1)=0.1476
c x(2)=0.5904
c sc=0.5
c x(1)=0.492
c x(2)=0.246
x(3)=0.195
x(4)=1.-x(1)-x(2)-x(3)
x(5)=1.
x(6)=0.
ff(7)=900+273
ff(8)=25+273
xco=0.
c write(11,9011)z,x(1),x(2),x(3),x(4)
c write(12,9011)z,x(5),x(6)
c write(13,9011)z,ff(7)-273,ff(8)-273
c write(14,9011)z,xco
c
c tube side
c
x10=x(1)
vv0t=vv1t*(ff(7)/298)*(1.01e5/ptt)
ft0t=ptt*vv0t/(rr*ff(7))
ut=vv0t/at
tt=at*ll/vv0t
c write(6,*) tt
c
c shell side
c
vv0s=vv1s*(ff(8)/298)*(1.01e5/pts)
ft0s=pts*vv0s/(rr*ff(8))
us=vv0s/as
ts=as*ll/vv0s
c write(6,9011) ft0t,ut,tt
keq=exp(4577.8/ff(7)-4.33)
cp(1)=(1.1-0.46*(ff(7)/1e3)+1.*(ff(7)/1e3)**2-
+0.454*(ff(7)/1e3)**3)*28.
cp(2)=(1.79-0.107*(ff(7)/1e3)+0.586*(ff(7)/1e3)**2-
+0.20*(ff(7)/1e3)**3)*18.
cp(3)=(13.46-4.6*(ff(7)/1e3)-6.85*(ff(7)/1e3)**2+
+3.79*(ff(7)/1e3)**3)*2.
cp(4)=(0.45+1.67*(ff(7)/1e3)-1.27*(ff(7)/1e3)**2+
+0.39*(ff(7)/1e3)**3)*44.
cp(5)=(5.193)*4.003
cp(6)=(13.46+4.6*(ff(8)/1e3)-6.85*(ff(8)/1e3)**2+
+3.79*(ff(8)/1e3)**3)*2.
rh=-4.12e4+(cp(3)+cp(4)-cp(1)-cp(2))*(ff(7)-298.)
do 201 i=1,4
ff(i)=ft0t*x(i)
ppt(i)=ptt*x(i)
201 continue
fcpt=0.
do 2011 i=1,4
fcpt=fcpt+ff(i)*cp(i)
2011 continue
do 202 i=5,6
ff(i)=ft0s*x(i)
ppt(i)=pts*x(i)
202 continue
fcps=0
do 2022 i=5,6
fcps=fcps+ff(i)*cp(i)
2022 continue
c L-H model
c
m=0
n=8
h=1e-5
zmax=11
1 beta=ptt(3)*ptt(4)/(ppt(1)*ppt(2)*keq)
rco=50*0.002093*
+(ppt(1)/(rr*ff(7)))**0.5*(ptt(2)/(rr*ff(7)))*(1.-beta)
pee=ppt(3)**0.5-ppt(6)**0.5
jh2=pmh2*pee
c jh2=0
keq=exp(4577.8/ff(7)-4.33)
cp(1)=(1.1-0.46*(ff(7)/1e3)+1.*(ff(7)/1e3)**2-
+0.454*(ff(7)/1e3)**3)*28.
cp(2)=(1.79-0.107*(ff(7)/1e3)+0.586*(ff(7)/1e3)**2-
+0.20*(ff(7)/1e3)**3)*18.
cp(3)=(13.46-4.6*(ff(7)/1e3)-6.85*(ff(7)/1e3)**2+
+3.79*(ff(7)/1e3)**3)*2.
cp(4)=(0.45+1.67*(ff(7)/1e3)-1.27*(ff(7)/1e3)**2+
+0.39*(ff(7)/1e3)**3)*44.
cp(5)=(5.193)*4.003
cp(6)=(13.46+4.6*(ff(8)/1e3)-6.85*(ff(8)/1e3)**2+
+3.79*(ff(8)/1e3)**3)*2.
rh=-4.12e4+(cp(3)+cp(4)-cp(1)-cp(2))*(ff(7)-298.)
f(1)=-at*rco
f(2)=-at*rco
f(3)=at*rco-jh2*pi*dt
f(4)=at*rco
f(5)=0.
f(6)=jh2*pi*dt
f(7)=(-rco*rh*at-pi*dt*jh2*cp(3)*(ff(7)-ff(8))-
+pi*dt*hts*(ff(7)-ff(8)))/fcpt
f(8)=(pi*dt*jh2*cp(3)*(ff(7)-298.)+pi*dt*hts*(ff(7)-ff(8))-
+pi*ds*hinf*(ff(8)-298.))/fcps
c f(8)=(pi*dt*hts*(ff(7)-ff(8))-pi*ds*hinf*(ff(8)-298))/fcps
call rks4(n,ff,f,z,h,m)
if(m.eq.5) go to 34
go to 1
c
c tube side
c
34 ft0t=0.
fcpt=0
do 101 i=1,4
ft0t=ft0t+ff(i)
fcpt=fcpt+ff(i)*cp(i)
101 continue
do 102 i=1,4
x(i)=ff(i)/ft0t
ppt(i)=ptt*x(i)
102 continue
c
c shell side
c
ft0s=0.
fcps=0
do 1011 i=5,6
ft0s=ft0s+ff(i)
fcps=fcps+ff(i)*cp(i)
1011 continue
do 1022 i=5,6
x(i)=ff(i)/ft0s
ppt(i)=pts*x(i)
1022 continue
if((z-zmax).ge.1e-6) go to 5
m=0
go to 1
5 xco=(x10-x(1))/x10*100
write(11, 9011) ptt/1.01e5,xco
if (ptt.gt.20*pref) go to 2002
ptt=ptt+5000.
go to 2001
2002 stop
9011 format(5e12.4)
1012 format(e12.4)
end
subroutine rks4(n,y,f,x,h,m)
implicit real*8(a-h,o-z)
dimemsion sum(100),savey(100),y(100),f(100)
m=m+1
go to (1,2,3,4,5),m
1 return
2 do 6 j=1,n
savey(j)=y(j)
sum(j)=f(j)
6 y(j)=savey(j)+0.5*h*f(j)
x=x+0.5*h
return
3 do 7 j=1,n
sum(j)=sum(j)+2.*f(j)
7 y(j)=savey(j)+0.5*h*f(j)
return
4 do 8 j=1,n
sum(j)=sum(j)+2.*f(j)
8 y(j)=savey(j)+h*f(j)
x=x+0.5*h
return
5 do 9 j=1,n
9 y(j)savey(j)+h*(sum(j)+f(j))/6.
return
end
有高手能幫忙嗎
報酬可議
--
※ 發信站: 批踢踢實業坊(ptt.cc)
◆ From: 125.230.82.61
1F：→ ripeSelf:codejob 版。 11/29 19:30
2F：→ geeoh:感謝 不知道有那個板 11/29 19:32
※ 發信站: 批踢踢實業坊(ptt.cc)
※ 轉錄者: geeoh (125.230.82.61), 時間: 11/29/2011 19:33:18
5F：→ fjm31714:好像有隻程是叫f2c 11/30 03:21
6F：→ WPC001:短短兩百行? 倒不如講一下你程式要做甚麼... 11/30 08:14
7F：→ WPC001:硬翻的話... 執行效率會不如fortran,理解的話,性能還OK 11/30 08:15
8F：→ WPC001:完全理解的話, C++也是有機會寫出比Fortran快的code 11/30 08:15
9F：推 sjgau:倒數 第三行的程式碼，似乎有錯？ 11/30 08:50