2010-11-04 98 views
0

我是一个非常糟糕的程序员,并且我被授予了一个程序来帮助我在空气动力学方面的工作。但它在Fortran中,而且我试图用MATLAB来运行这个程序。任何帮助将其转换为语言matlab理解? (preferabbly C++)需要将以下FORTRAN代码转换为C++

 program joukow 
c 

c computes joukowski airfoil and finds pressure coefficient 

c currently set up for symmetric airfoil with sharp trailing edge 

c and chord length equal to one. 

c profile is written onto prof.dat and cp onto cp.dat 

c  implicit real*8(a-h,o-z) 

     complex z,zeta,cw 
     dimension uz(100),vz(100),xi(100),eta(100),cp(100) 
     dimension xout(100),yout(100) 
     open(unit=8,file='prof.dat',status='unknown') 
     open(unit=9,file='cp.dat',status='unknown') 
     b=1.d0 
     write(6,98) 
     format(2x,'input the radius of the a-circle in z plane') 
     read(5,99)a 
     format(f10.0) 
     xl=2.*a-1.+1./(2.*a-1.) 

c  xl=a+1./a 

c  chord=2.*xl 

     chord=2.+xl 
     del=a-b 

c  del =0.1d0 
     do 50 i=1,100 
     ri=i 
     theta=6.2832d0*ri/101.d0 
     x=-del+a*cos(theta) 
     y=a*sin(theta) 
     z=cmplx(x,y) 
     zeta=z+b**2/z 

c 

c xi and eta are coordinates of points on airfoil 

c 

     xi(i)=real(zeta) 

     eta(i)=aimag(zeta) 

     cw=(1.-a**2/(z+del)**2)/(1.-b**2/z**2) 
c 

c uz and vz are velocity components on the airfoil assuming the free-stream 

c speed is one. 
c 
     uz(i)=real(cw) 
     vz(i)=-aimag(cw) 

c 

c xout and yout are airfoil coordinates where the leading edge is at (0,0) 

c and the chordlength is one. 

c 

     xout(i)=(xl+xi(i))/chord 
     yout(i)=eta(i)/chord 
     write(8,100)xout(i),yout(i) 
     format(2x,2f10.4) 
     continue 

c 

c now calculate the pressure coefficient cp 

c 

     write(6,200) 
     format(2x,'pressure coefficients') 
     do 70 i=1,50 
     cp(i)=1.-(uz(i)**2+vz(i)**2) 
     write(9,100)xout(i),cp(i) 
     continue 
     stop 
     end 
+1

您提交的程序似乎是不完整的 - 没有像'do 70 i = 1,50'这样的语句编号(参考语句编号70作为循环的结尾)。 – 2010-11-04 17:06:35

回答

0

我不知道它是如何以及仍然支持 - - 而是直接转化成FORTRAN C代码曾经是F2C最简单的方法。

+0

我不知道你是否尝试过在F77和更老的版本上使用它,但它有点痛苦。上次我尝试使用它时,我不得不完全重写FORTRAN。 – greyfade 2010-11-04 17:25:50

+0

我不得不同意greyface,f2c作为编译过程的一部分是很好的,但输出不适合人类消费。 – dmckee 2010-11-04 17:36:56

+0

但是,如果你只是把Fortran转换成一个matlab库,你关心它的可读性如何? – 2010-11-04 17:48:57

8

Matlab理解Fortran就好 - 检查文档。如果这不能满足您的需求,程序中的大部分线路都可以通过很少的修改就可以输入Matlab控制台。如果你是一个糟糕的程序员,我建议你花时间将程序修改为Matlab而不是C++。如果你没有得到比我现在有空的更好的帮助,我会写更多的。

编辑:首先,关于来自Matlab的using Fortran source files的一些信息。如果你真的不想(或者不能或者没有性能原因,不要这样做)将Fortran重写为Matlab,然后将其转换为MEX文件。使用f2c(或其他任何东西,包括您自己的时间和精力)首先将Fortran转换为C或C++对我来说似乎毫无意义。

如果你不喜欢这个想法,下面是将Fortran转换为Matlab的一些想法。首先,所有以C或c开头的行都是注释,所以你不需要翻译它们。从您的代码开始:

complex z,zeta,cw 
    dimension uz(100),vz(100),xi(100),eta(100),cp(100) 
    dimension xout(100),yout(100) 

这些行声明了一些变量。在Matlab中使用它们之前,您不必声明变量,但有时候有很好的理由这样做。你也不必在Fortran中,尽管这些日子普遍被认为是一个坏主意。你可以“声明”在Matlab这些变量与语句如:

uz = zeros(100,1); 
vz = zeros(100,1); 

通过预先在您的MATLAB声明这些你为他们一次分配内存,并避免一些性能降低的问题。

接下来的2行:

 open(unit=8,file='prof.dat',status='unknown') 
    open(unit=9,file='cp.dat',status='unknown') 

打开输出文件的情侣。它们稍后在write语句中使用 - 忘记它们,而是编写Matlab语句,例如save xout

下一行是Fortran语言,但同样在Matlab:

b=1.d0 

下一行获取值从控制台半径:

write(6,98) 
    format(2x,'input the radius of the a-circle in z plane') 
    read(5,99)a 
    format(f10.0) 

再次,我建议你忘记这些,只是使用Matlab控制台设置值为a。更多的Fortran不需要被翻译(尽管我建议你或者在没有跟随0的情况下删除小数点,或者在它们之间放置一个空格并且随后的* - 。*是Matlab中的特定操作符):

xl=2.*a-1.+1./(2.*a-1.) 

    chord=2.+xl 
    del=a-b 

Fortran do循环与Matlab for循环相同。重写:

do 50 i=1,100 

for i = 1:100 

由于其他受访者指出,其中相匹配的结束语句超出目前还不清楚,你必须明白这一点。请注意,我只是提供了Fortran到Matlab的逐行转换。它不是精心编写的Fortran,而且我没有提供写得很好的Matlab,我会把它留给你。

这很多并不需要翻译:

ri=i 
    theta=6.2832d0*ri/101.d0 
    x=-del+a*cos(theta) 
    y=a*sin(theta) 

CMPLX是Fortran函数返回一个复数具有实部x和虚部Y:

z=cmplx(x,y) 

在Matlab中这将是z = x + y * i。 Fortran使用**进行取幂,Matlab使用^

zeta=z+b**2/z 

等等。

希望有所帮助。

2

我用f2matlab,稍后摸了一下。这里是清理和编译代码FORTRAN90:

program joukow 
! 
! computes joukowski airfoil and finds pressure coefficient 
! currently set up for symmetric airfoil with sharp trailing edge 
! and chord length equal to one. 
! profile is written onto prof.dat and cp onto cp.dat 
!  implicit real*8(a-h,o-z) 
complex z,zeta,cw 
dimension uz(100),vz(100),xi(100),eta(100),cp(100) 
dimension xout(100),yout(100) 
open(unit=8,file='prof.dat',status='unknown') 
open(unit=9,file='cp.dat',status='unknown') 
b=1.d0 
write(6,98) 
98 format(2x,'input the radius of the a-circle in z plane') 
read(5,99)a 
99 format(f10.0) 
xl=2.*a-1.+1./(2.*a-1.) 
!  xl=a+1./a 
!  chord=2.*xl 
chord=2.+xl 
del=a-b 
!  del =0.1d0 
do i=1,100 
    ri=i 
    theta=6.2832d0*ri/101.d0 
    x=-del+a*cos(theta) 
    y=a*sin(theta) 
    z=cmplx(x,y) 
    zeta=z+b**2/z 
    ! 
    ! xi and eta are coordinates of points on airfoil 
    ! 
    xi(i)=real(zeta) 
    eta(i)=aimag(zeta) 
    cw=(1.-a**2/(z+del)**2)/(1.-b**2/z**2) 
    ! 
    ! uz and vz are velocity components on the airfoil assuming the free-stream 
    ! speed is one. 
    ! 
    uz(i)=real(cw) 
    vz(i)=-aimag(cw) 
    ! 
    ! xout and yout are airfoil coordinates where the leading edge is at (0,0) 
    ! and the chordlength is one. 
    ! 
    xout(i)=(xl+xi(i))/chord 
    yout(i)=eta(i)/chord 
    write(8,100)xout(i),yout(i) 
100 format(2x,2f10.4) 
end do 
! 
! now calculate the pressure coefficient cp 
! 
write(6,200) 
200 format(2x,'pressure coefficients') 
do i=1,50 
    cp(i)=1.-(uz(i)**2+vz(i)**2) 
    write(9,100) xout(i),cp(i) 
end do 
stop 
end program joukow 

这里是产生MATLAB代码:

function hw1(varargin) 
% 
% computes joukowski airfoil and finds pressure coefficient 
% currently set up for symmetric airfoil with sharp trailing edge 
% and chord length equal to one. 
% profile is written onto prof.dat and cp onto cp.dat 
%  implicit real*8(a-h,o-z) 

format_99=['%10.0f']; 
format_100=[repmat(' ',1,2),repmat('%10.4f',1,2),'\n']; 
format_200=[repmat(' ',1,2),'pressure coefficients \n']; 

fid_8=fopen('prof.dat','w+'); 
fid_9=fopen('cp.dat','w+'); 
b=1.0d0; 
a=input('input the radius of the a-circle in z plane'); 
xl=2..*a-1.+1../(2..*a-1.); 
%  xl=a+1./a 
%  chord=2.*xl 
chord=2.+xl; 
del=a-b; 
%  del =0.1d0 
for i=1:100; 
ri=i; 
theta=6.2832d0.*ri./101.0d0; 
x=-del+a.*cos(theta); 
y=a.*sin(theta); 
z=complex(x,y); 
zeta=z+b.^2./z; 
% 
% xi and eta are coordinates of points on airfoil 
% 
xi(i)=real(zeta); 
eta(i)=imag(zeta); 
cw=(1.-a.^2./(z+del).^2)./(1.-b.^2./z.^2); 
% 
% uz and vz are velocity components on the airfoil assuming the free-stream 
% speed is one. 
% 
uz(i)=real(cw); 
vz(i)=-imag(cw); 
% 
% xout and yout are airfoil coordinates where the leading edge is at (0,0) 
% and the chordlength is one. 
% 
xout(i)=(xl+xi(i))./chord; 
yout(i)=eta(i)./chord; 
fprintf(fid_8,format_100,xout(i),yout(i)); 
end; i=100+1; 
% 
% now calculate the pressure coefficient cp 
% 
fprintf(1,format_200); 
for i=1:50; 
cp(i)=1.-(uz(i).^2+vz(i).^2); 
fprintf(fid_9,format_100, xout(i),cp(i)); 
end; i=50+1; 
end %program joukow 

他们都给出了相同的结果对我来说。但是,我没有检查算法的正确性,只是转换了代码。