2010-11-04 108 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 

他們都給出了相同的結果對我來說。但是,我沒有檢查算法的正確性,只是轉換了代碼。