Необходимо преобразовать следующий код FORTRAN в C ++

Question

Я очень плохой программист, и мне дали программу, которая якобы поможет мне в моей аэродинамике. но это в Фортране, и я пытаюсь использовать MATLAB для запуска этой программы. любая помощь по преобразованию его в язык, который понимает Matlab? (предпочтительно с ++)

      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

Answer 1

Matlab прекрасно понимает Фортран - проверьте документацию.И если вас это не устраивает, большинство строк в программе, которые выполняют какие-либо вычисления, могут быть напечатаны в консоли Matlab с очень небольшими изменениями.Если вы плохой программист, я советую вам потратить время на модификацию программы в Matlab, а не в C ++.Я напишу позже, если вы не получите никакой лучшей помощи, чем у меня есть время.

РЕДАКТИРОВАТЬ: во-первых, некоторая информация о с использованием исходных файлов Fortran от Matlab.Если вы действительно не хотите (или не можете, или у вас есть причины не делать этого) переписать Фортран в 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, но иногда есть веские причины для этого.Вам не нужно делать это на Фортране, хотя в наши дни это повсеместно считается плохой идеей.Вы можете «объявить» эти переменные в 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.Больше Фортрана, который не нужно переводить (хотя я предлагаю вам либо отбросить десятичные точки без следующих 0, либо поставить пробел между ними и последующим * -. * - это специальный оператор в Matlab):

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

  chord=2.+xl
  del=a-b

Цикл do Фортрана такой же, как цикл Matlab for.Перепишите:

  do 50 i=1,100

как

for i = 1:100

Как заметил один из респондентов, неясно, куда идет соответствующее конечное утверждение, вам придется это выяснить.Обратите внимание, что я просто предлагаю построчный перевод Фортрана на Matlab.Это не хорошо написанный Фортран, и я не предлагаю хорошо написанный Matlab, я оставлю это вам.

Этот лот не нужно переводить:

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

cmplx - это функция Фортрана, которая возвращает комплексное число, которое имеет вещественную часть x и мнимую часть y:

  z=cmplx(x,y)

В Matlab это будет z = x + y * i.Fortran использует ** для возведения в степень, Matlab использует ^

  zeta=z+b**2/z

и т. Д. И т. П.

Надеюсь, что это поможет.

Answer 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

Они оба дают одинаковые результаты для меня.Я не проверял алгоритм на корректность, просто преобразовал код.

Answer 3

Я не знаю, насколько хорошо он все еще поддерживается - но самым простым способом был f2c, который переводит фортран непосредственно в код на языке c.