help to convert this fortan program to c++

Discussion in 'C++' started by varunkumarchundayil@gmail.com, Apr 11, 2009.

  1. Guest

    C PROGRAMME TROJANS
    PROGRAM NBODY
    C THIS IS A GENERAL N-BODY PROGRAMME WHERE INTER-BODY FORCES BETWEEN
    C BODIES i AND j ARE OF THE FORM CM(i)*CM(j)*F(Rij,Vij) WHERE F IS A
    C FUNCTION OF Rij, THE DISTANCE BETWEEN THE BODIES, AND Vij, THE
    C RELATIVE VELOCITIES OF THE TWO BODIES.
    C THE STRUCTURE OF THE PROGRAMME IS;
    C
    C (I) THE MAIN PROGRAMME "NBODY"IS WHICH INCLUDES THE RUNGE-KUTTA
    C ROUTINE WITH AUTOMATIC STEP CONTROL.
    C (II) SUBROUTINE "START" WHICH ENABLES INPUT OF THE INITIAL
    BOUNDARY
    C CONDITIONS.
    C (III) SUBROUTINE "ACC" WHICH GIVES THE ACCELERATION OF EACH BODY
    C DUE TO ITS INTERACTIONS WITH ALL OTHER BODIES.
    C (IV) SUBROUTINE "STORE" WHICH STORES INTERMEDIATE COORDINATES AND
    C VELOCITY COMPONENTS AS THE COMPUTATION PROGRESSES.
    C (V) SUBROUTINE "OUT" WHICH OUTPUTS THE RESULTS TO DATA FILES.
    C
    C BY CHANGING THE SUBROUTINES DIFFERENT PROBLEMS MAY BE SOLVED.
    C THE CM'S CAN BE MASSES OR CHARGES OR BE MADE EQUAL TO UNITY WHILE
    C THE FORCE LAW CAN BE INVERSE-SQUARE OR ANYTHING ELSE -
    C e.g. LENNARD-JONES. SEE COMMENT AT THE BEGINNING OF SUBROUTINE
    C "ACC" FOR THE TYPES OF FORCES OPERATING.
    C
    C THE FOUR-STEP RUNGE-KUTTA ALGORITHM IS USED. THE RESULTS OF TWO
    C STEPS WITH TIMESTEP H ARE CHECKED AGAINST TAKING ONE STEP WITH
    C TIMESTEP 2*H. IF THE DIFFERENCE IS WITHIN THE TOLERANCE THEN THE
    C TWO STEPS, EACH OF H, ARE ACCEPTED AND THE STEPLENGTH IS DOUBLED
    FOR
    C THE NEXT STEP. HOWEVER, IF THE TOLERANCE IS NOT SATISFIED THEN
    THE
    C STEP IS NOT ACCEPTED AND ONE TRIES AGAIN WITH A HALVED STEPLENGTH.
    C IT IS ADVISABLE, BUT NOT ESSENTIAL, TO START WITH A REASONABLE
    C STEPLENGTH; THE PROGRAMME QUICKLY FINDS A SUITABLE VALUE.
    C
    C AS PROVIDED THE PROGRAMME HANDLES UP TO 20 BODIES BUT THIS CAN BE
    C CHANGED FROM 20 TO WHATEVER IS REQUIRED IN THE DIMENSION STATEMENT.
    C
    DIMENSION CM(20),X(20,3),V(20,3),DX(20,3,0:4),DV(20,3,0:4),WT
    (4),
    +XTEMP(2,20,3),VTEMP(2,20,3),XT(20,3),VT(20,3),DELV(20,3)
    COMMON/A/X,V,TOL,H,TOTIME,DELV,XT,VT,NB,IST,TIME,IG,XTEMP,
    +VTEMP,CM
    DATA WT/0.0,0.5,0.5,1.0/
    IST=0
    OPEN(UNIT=9,FILE='LPT1')
    C
    C SETTING THE INITIAL BOUNDARY CONDITION CAN BE DONE EITHER
    C EXPLICITLY AS VALUES OF (X,Y,Z) AND (U,V,W) FOR EACH BODY OR
    C CAN BE COMPUTED. THIS IS CONTROLLED BY SUBROUTINE "START".
    C OTHER PARAMETERS ARE ALSO SET IN "START" WHICH ALSO INDICATES
    C THE SYSTEM OF UNITS BEING USED.
    C
    CALL START
    C
    TIME=0
    C INITIALIZE ARRAYS
    DO 57 I=1,20
    DO 57 J=1,3
    DO 57 K=0,4
    DX(I,J,K)=0
    DV(I,J,K)=0
    57 CONTINUE
    C
    C WE NOW TAKE TWO STEPS WITH STEP LENGTH H FOLLOWED BY ONE STEP
    C WITH STEP LENGTH 2*H FROM THE SAME STARTING POINT BUT FIRST WE
    C STORE THE ORGINAL SPACE AND VELOCITY COORDINATES AND TIMESTEP.
    C
    25 DO 7 IT=1,2
    DO 8 J=1,NB
    DO 8 K=1,3
    XTEMP(IT,J,K)=X(J,K)
    VTEMP(IT,J,K)=V(J,K)
    8 CONTINUE
    HTEMP=H
    DO 10 NOSTEP=1,3-IT
    DO 11 I=1,4
    DO 12 J=1,NB
    DO 12 K=1,3
    XT(J,K)=XTEMP(IT,J,K)+WT(I)*DX(J,K,I-1)
    12 VT(J,K)=VTEMP(IT,J,K)+WT(I)*DV(J,K,I-1)
    C
    CALL ACC
    C
    DO 13 J=1,NB
    DO 13 K=1,3
    DV(J,K,I)=IT*HTEMP*DELV(J,K)
    DX(J,K,I)=IT*HTEMP*VT(J,K)
    13 CONTINUE
    11 CONTINUE
    DO 14 J=1,NB
    DO 14 K=1,3
    XTEMP(IT,J,K)=XTEMP(IT,J,K)+(DX(J,K,1)+DX(J,K,4)+2*
    +(DX(J,K,2)+DX(J,K,3)))/6.0
    VTEMP(IT,J,K)=VTEMP(IT,J,K)+(DV(J,K,1)+DV(J,K,4)+2*
    +(DV(J,K,2)+DV(J,K,3)))/6.0
    14 CONTINUE
    10 CONTINUE
    7 CONTINUE
    C
    C THE ABOVE HAS MADE TWO STEPS OF H AND, FROM THE SAME STARTING
    POINT,
    C A SINGLE STEP OF 2*H. THE RESULTS ARE NOW COMPARED
    C
    DO 20 J=1,NB
    DO 20 K=1,3
    IF(ABS(XTEMP(1,J,K)-XTEMP(2,J,K)).GT.TOL)THEN
    H=0.5*H
    GOTO 25
    ENDIF
    20 CONTINUE
    C
    C AT THIS STAGE THE DOUBLE STEP WITH H AGREES WITHIN TOLERANCE WITH
    C THE SINGLE STEP WITH 2*H. THE TIMESTEP WILL NOW BE TRIED WITH
    C TWICE THE VALUE FOR THE NEXT STEP. IF IT IS TOO BIG THEN IT WILL
    C BE REDUCED AGAIN.
    C
    H=2*H
    DO 80 J=1,NB
    DO 80 K=1,3
    X(J,K)=XTEMP(1,J,K)
    V(J,K)=VTEMP(1,J,K)
    80 CONTINUE
    TIME=TIME+H
    C
    CALL STORE
    C
    IF(TIME.GE.TOTIME)GOTO 50
    IF(IG.GT.1000)THEN
    IG=1000
    GOTO 50
    ENDIF
    GOTO 25
    C
    50 CALL OUT
    C
    STOP
    END


    SUBROUTINE STORE
    DIMENSION CM(20),X(20,3),V(20,3),XSTORE(1000,20,3),
    +XTEMP(2,20,3),VTEMP(2,20,3),DELV(20,3),XT(20,3),VT(20,3)
    COMMON/A/X,V,TOL,H,TOTIME,DELV,XT,VT,NB,IST,TIME,IG,XTEMP,
    +VTEMP,CM
    COMMON/B/NORIG,XSTORE
    DO 21 J=1,NB
    DO 21 K=1,3
    X(J,K)=XTEMP(1,J,K)
    V(J,K)=VTEMP(1,J,K)
    21 CONTINUE
    C
    C UP TO 1000 POSITIONS ARE STORED. THESE ARE TAKEN EVERY 50 STEPS.
    C
    IST=IST+1
    IF((IST/50)*50.NE.IST)GOTO 50
    IG=IST/50
    IF(IG.GT.1000)GOTO 50
    DO 22 J=1,NB
    DO 22 K=1,3
    XSTORE(IG,J,K)=X(J,K)
    22 CONTINUE
    50 RETURN
    END


    SUBROUTINE START
    DIMENSION CM(20),X(20,3),V(20,3),XSTORE(1000,20,3),
    +XTEMP(2,20,3),VTEMP(2,20,3),DELV(20,3),XT(20,3),VT(20,3)
    COMMON/A/X,V,TOL,H,TOTIME,DELV,XT,VT,NB,IST,TIME,IG,XTEMP,
    +VTEMP,CM
    COMMON/B/NORIG,XSTORE
    OPEN(UNIT=21,FILE='LAST.DAT')
    OPEN(UNIT=22,FILE='FAST.DAT')
    C
    C THE PROGRAMME AS PROVIDED IS FOR AN INVERSE-SQUARE LAW AND
    C USES UNITS FOR WHICH THE UNIT MASS IS THAT OF THE SUN, THE
    C THE UNIT OF DISTANCE IS THE ASTRONOMICAL UNIT (MEAN SUN-EARTH
    C DISTANCE),THE UNIT OF TIME IS THE YEAR AND THE GRAVITATIONAL
    C CONSTANT IS 4*PI**2.
    C
    WRITE(6,'('' INPUT THE NUMBER OF BODIES'')')
    READ(5,*)NB
    WRITE(6,'('' INPUT THE VALUES OF CM IN SOLAR-MASS UNITS. '')')
    WRITE(6,'('' THE TROJAN ASTEROID MASSES CAN BE PUT AS ZERO'')')
    DO 1 I=1,NB
    WRITE(6,500)I
    500 FORMAT(25H READ IN THE VALUE OF CM[,I3, 1H])
    READ(5,*)CM(I)
    1 CONTINUE
    WRITE(6,'('' INPUT THE INITIAL TIMESTEP [years]'')')
    READ(5,*)H
    WRITE(6,'('' INPUT TOTAL TIME FOR THE SIMULATION [years]'')')
    READ(5,*)TOTIME
    C
    C THE PROGRAMME ASKS THE USER TO SPECIFY A TOLERANCE, THE MAXIMUM
    C ABSOLUTE ERROR THAT CAN BE TOLERATED IN ANY POSITIONAL COORDINATE
    C (X, Y OR Z). IF THIS IS SET TOO LOW THEN THE PROGRAMME CAN BECOME
    C VERY SLOW. FOR COMPUTATIONS INVOLVING PLANETS A TOLERANCE OF
    1.0E-6
    C (c. 150 KM) IS USUALLY SATISFACTORY.
    C
    WRITE(6,'('' INPUT THE TOLERANCE '')')
    WRITE(6,'('' SEE COMMENT ABOVE THIS STATEMENT IN LISTING'')')
    READ(5,*)TOL
    WRITE(6,'('' THE CALCULATION CAN BE DONE RELATIVE TO AN '')')
    WRITE(6,'('' ARBITRARY ORIGIN OR WITH RESPECT TO ONE OF '')')
    WRITE(6,'('' THE BODIES AS ORIGIN. INPUT ZERO FOR AN '')')
    WRITE(6,'('' ARBITRARY ORIGIN OR THE NUMBER OF THE BODY.'')')
    WRITE(6,'('' IF A BODY IS CHOSEN AS ORIGIN THEN ALL ITS'')')
    WRITE(6,'('' POSITIONAL AND VELOCITY VALUES ARE SET TO
    ZERO'')')
    READ(5,*)NORIG
    DO 31 J=1,NB
    WRITE(6,100)J
    100 FORMAT(23H INPUT [X,Y,Z] FOR BODY,I3)
    READ(5,*)X(J,1),X(J,2),X(J,3)
    WRITE(6,200)J
    200 FORMAT(32H INPUT [XDOT,YDOT,ZDOT] FOR BODY,I3)
    READ(5,*)V(J,1),V(J,2),V(J,3)
    31 CONTINUE
    RETURN
    END




    SUBROUTINE ACC
    DIMENSION CM(20),X(20,3),V(20,3),XSTORE(1000,20,3),R(3),
    +XTEMP(2,20,3),VTEMP(2,20,3),DELV(20,3),XT(20,3),VT(20,3),DD(3)
    COMMON/A/X,V,TOL,H,TOTIME,DELV,XT,VT,NB,IST,TIME,IG,XTEMP,
    +VTEMP,CM
    COMMON/B/NORIG,XSTORE
    C
    C THE PROGRAMME AS PROVIDED IS FOR AN INVERSE-
    C SQUARE LAW AND USES UNITS FOR WHICH THE UNIT MASS IS THAT OF THE
    SUN,
    C THE UNIT OF DISTANCE IS THE ASTRONOMICAL UNIT (MEAN SUN-EARTH
    DISTANCE),
    C THE UNIT OF TIME IS THE YEAR AND THE GRAVITATIONAL CONSTANT IS
    4*PI**2.
    C HOWEVER, THE USER MAY MODIFY THE SUBROUTINE "ACC" TO CHANGE TO ANY
    OTHER
    C FORCE LAW AND/OR ANY OTHER SYSTEM OF UNITS.
    C
    C SET THE VALUE OF G IN ASTRONOMICAL UNITS
    PI=4.0*ATAN(1.0)
    G=4*PI*PI
    DO 1 J=1,NB
    DO 1 K=1,3
    DELV(J,K)=0
    1 CONTINUE
    C THE FOLLOWING PAIR OF DO LOOPS FINDS INTERACTIONS FOR ALL PAIRS
    C OF BODIES
    DO 2 J=1,NB-1
    DO 2 L=J+1,NB
    DO 3 K=1,3
    R(K)=XT(J,K)-XT(L,K)
    3 CONTINUE
    RRR=(R(1)**2+R(2)**2+R(3)**2)**1.5
    DO 4 K=1,3
    C THE NEXT TWO STATEMENTS GIVE THE CONTRIBUTIONS TO THE THREE
    C COMPONENTS OF ACCELERATION DUE TO BODY J ON BODY L AND DUE TO
    C BODY L ON BODY J.
    DELV(J,K)=DELV(J,K)-G*CM(L)*R(K)/RRR
    DELV(L,K)=DELV(L,K)+G*CM(J)*R(K)/RRR
    4 CONTINUE
    2 CONTINUE
    C IF ONE OF THE BODIES IS TO BE THE ORIGIN THEN ITS ACCELERATION IS
    C SUBTRACTED FROM THAT OF ALL OTHER BODIES.
    IF(NORIG.EQ.0)GOTO 10
    DD(1)=DELV(NORIG,1)
    DD(2)=DELV(NORIG,2)
    DD(3)=DELV(NORIG,3)
    DO 6 J=1,NB
    DO 6 K=1,3
    DELV(J,K)=DELV(J,K)-DD(K)
    6 CONTINUE
    10 RETURN
    END


    SUBROUTINE OUT
    DIMENSION CM(20),X(20,3),V(20,3),XSTORE(1000,20,3),
    +XTEMP(2,20,3),VTEMP(2,20,3),DELV(20,3),XT(20,3),VT(20,3)
    COMMON/A/X,V,TOL,H,TOTIME,DELV,XT,VT,NB,IST,TIME,IG,XTEMP,
    +VTEMP,CM
    COMMON/B/NORIG,XSTORE
    C (X,Y) VALUES FOR THE LEADING ASTEROID ARE PLACED IN DATA FILE
    C LAST.DAT AND FOR THE FOLLOWING ASTEROID IN DATA FILE FAST.DAT.
    C
    C FOR THE TROJAN ASTEROID PROBLEM THE JUPITER RADIUS VECTOR IS
    C ROTATED TO PUT IT ON THE Y AXIS. THE POSITIONS OF THE ASTEROIDS
    C RELATIVE TO JUPITER ARE PLOTTED.
    LAST=MIN(IG,1000)
    DO 60 I=1,LAST
    C THETA IS THE ANGLE BETWEEN THE X AXIS AND THE JUPITER RADIUS VECTOR
    THETA=ATAN2(XSTORE(I,2,2),XSTORE(I,2,1))
    XSTORE(I,2,2)=SQRT(XSTORE(I,2,2)**2+XSTORE(I,2,1)**2)
    XSTORE(I,2,1)=0
    C NOW THE ASTEROID RADIUS VECTORS ARE ROTATED BY PI/2-THETA
    DO 61 J=1,2
    AA=XSTORE(I,J+2,1)*SIN(THETA)-XSTORE(I,J+2,2)*COS(THETA)
    BB=XSTORE(I,J+2,1)*COS(THETA)+XSTORE(I,J+2,2)*SIN(THETA)
    XSTORE(I,J+2,1)=AA
    XSTORE(I,J+2,2)=BB
    61 CONTINUE
    60 CONTINUE
    C THE MODIFIED POSITIONS ARE NOW OUTPUT TO DATA FILES.
    DO 63 J=3,4
    N=18+J
    REWIND N
    DO 64 I=1,LAST
    WRITE(N,*)XSTORE(I,J,1),XSTORE(I,J,2)
    64 CONTINUE
    63 CONTINUE
    RETURN
    END
     
    , Apr 11, 2009
    #1
    1. Advertising

  2. red floyd Guest

    wrote:
    [redacted]

    How much are you planning on paying me to do so?
     
    red floyd, Apr 11, 2009
    #2
    1. Advertising

  3. osmium Guest

    <> wrote:


    >C PROGRAMME TROJANS
    > PROGRAM NBODY
    > C THIS IS A GENERAL N-BODY PROGRAMME WHERE INTER-BODY FORCES BETWEEN
    > C BODIES i AND j ARE OF THE FORM CM(i)*CM(j)*F(Rij,Vij) WHERE F IS A
    > C FUNCTION OF Rij, THE DISTANCE BETWEEN THE BODIES, AND Vij, THE
    > C RELATIVE VELOCITIES OF THE TWO BODIES.
    > C THE STRUCTURE OF THE PROGRAMME IS;

    <huge snip>

    That is pretty nasty looking stuff by today's standards. I suggest looking
    in _Numerical Recipes_ (in C). That code, or an equivalent, might be in
    there. I guess this is the kind of thing that mandated Engineering majors
    were required to study Fortran.

    Note that you asked for help but didn't provide a start. IOW you want to
    farm out the whole darn thing!
     
    osmium, Apr 11, 2009
    #3
    1. Advertising

Want to reply to this thread or ask your own question?

It takes just 2 minutes to sign up (and it's free!). Just click the sign up button to choose a username and then you can ask your own questions on the forum.
Similar Threads
  1. ola
    Replies:
    3
    Views:
    714
    Marco Schmidt
    Feb 16, 2004
  2. tomekj
    Replies:
    1
    Views:
    383
    stelios xanthakis
    Oct 29, 2003
  3. tomekj
    Replies:
    2
    Views:
    390
    tomekj
    Oct 29, 2003
  4. Coca
    Replies:
    15
    Views:
    650
    Alan Balmer
    Jan 14, 2004
  5. Alex van der Spek

    Calling FORTAN dll functions from Python

    Alex van der Spek, Dec 7, 2010, in forum: Python
    Replies:
    3
    Views:
    461
    Carl Banks
    Dec 8, 2010
Loading...

Share This Page