The re-engineered code is shown below. The changes are:
A header has been added to record when and how the code was processed;
A USE statement has been inserted for the emulation module;
The declarations of the real variables have been changed to the emulated derived type em_real_k4;
The parameter values have been converted to objects of the derived type em_real_k4;
!H!**************************************************************************** !H! File: ../fpt_output/ss_9b.f90 !H! Output by FPT 3.8-j Intel-Linux On 24-SEP-15 At 16:53:47 Input files: !H! Main: /media/truecrypt1/fpt/emulated_reals/ss/fpt/ss_9b_precision.fsp !H! Current: /media/truecrypt1/fpt/emulated_reals/ss/original_source/ss_9b.f90 !H! Licensee: SimCon: Development version. !H!**************************************************************************** ! ss_9b.f90 16-Sep-13 John Collins ! 16-Sep-13 John Collins ! PROGRAM ss_9_body ! ! USE module_emulate_real_arithmetic ! IMPLICIT NONE ! INTEGER,PARAMETER :: & kr = KIND(1.0D0) ! Real kind INTEGER,PARAMETER :: & frames_per_day = 32 ! INTEGER,PARAMETER :: & sample_interval = 30*frames_per_day ! 30 days INTEGER,PARAMETER :: & end_frame = (10*365+2)* & frames_per_day ! ! Note: We measure time in days, distance in AU, mass in kg ! TYPE(em_real_k8),PARAMETER :: & au = em_real_k8(149597870700.0D0) TYPE(em_real_k8),PARAMETER :: & gu = em_real_k8(-6.67300D-11) ! m**3/kg/s**2 TYPE(em_real_k8),PARAMETER :: & g = em_real_k8(((60.0D0*60.0D0*24.0D0)**2/au%value**3)*gu%value) ! au**3/kg/day**2 TYPE(em_real_k8) & time & ,p_time & ! Previous time ,dt & ,end_time & ,mass(9) & ,s(1:3,1:9) & ! Position (axis,body) ,ds(1:3,1:9,1:9) & ! Difference in position, (axis,body1,body2) ,s2(1:3) & ! Square of difference in position, body pair ,r2(1:9,1:9) & ! Square of radial distance, body to body ,r(1:9,1:9) & ! Radial distances, body to body ,gf(1:9,1:9) & ! Total gravitational field, body to body ,proj(1:3,1:9,1:9) & ! Cosine projections on axes (axis,body1,body2) ,vf(1:3,1:9,1:9) & ! Vector field components (axis,body1,body2) ,a(1:3,1:9) & ! Accelerations (axis,body) ,pv(1:3,1:9) & ! Previous velocities (axis,body) ,v(1:3,1:9) & ! Velocities (axis,body) ,pa(1:3,1:9) ! Previous accelerations (axis,body) ! INTEGER & axis & , body1 & , body2 & , frame ! LOGICAL & end_run ! DATA & mass /em_real_k8(1.988544D+30) & ! Sun ,em_real_k8(3.302D+23) & ! Mercury ,em_real_k8(48.685D+23) & ! Venus ,em_real_k8(59.7219D+23) & ! Earth ,em_real_k8(6.4185D+23) & ! Mars ,em_real_k8(18981.3D+23) & ! Jupiter ,em_real_k8(5683.19D+23) & ! Saturn ,em_real_k8(868.103D+23) & ! Uranus ,em_real_k8(1024.1D+23)/ & ! Neptune ! ! Initialisation data from JPL, California Institute of Technology, http://ssd.jpl.nasa.gov/horizons.cgi#top ! Positions in AU ,s /em_real_k8(-3.747147775505799D-03),em_real_k8(2.926587035303590D-03),em_real_k8(4.446807296978120D-06) & ! Sun ,em_real_k8(4.715831111175716D-02),em_real_k8(3.056015767171350D-01),em_real_k8(2.006272440161816D-02) & ! Mercury ,em_real_k8(4.959853421725741D-02),em_real_k8(-7.222408775972484D-01),em_real_k8(-1.300479358383397D-02) & ! Venus ,em_real_k8(-1.797649390463135D-01),em_real_k8(9.703470886653921D-01),em_real_k8(-1.743826995752289D-05) & ! Earth ,em_real_k8(-7.334184274376777D-01),em_real_k8(1.457212408029754D+00),em_real_k8(4.839432410544958D-02) & ! Mars ,em_real_k8(4.505316309281244D+00),em_real_k8(-2.163555523180046D+00),em_real_k8(-9.190499900147728D-02) & ! Jupiter ,em_real_k8(-9.468003745121990D+00),em_real_k8(2.616109582176896D-01),em_real_k8(3.722149420390660D-01) & ! Saturn ,em_real_k8(2.003287446208430D+01),em_real_k8(-1.529803159589830D+00),em_real_k8(-2.652256724942619D-01) & ! Uranus ,em_real_k8(2.481340360454546D+01),em_real_k8(-1.689465149553225D+01),em_real_k8(-2.239358924631250D-01)/ & ! Neptune ! ! Velocities in AU per day ,v/em_real_k8(-2.990258994506073D-06),em_real_k8(-5.530878190433255D-06),em_real_k8(6.981801439481427D-08) & ! Sun ,em_real_k8(-3.338934398882634D-02),em_real_k8(5.702520017360675D-03),em_real_k8(3.529925953990557D-03) & ! Mercury ,em_real_k8(2.003227503747324D-02),em_real_k8(1.405250683890825D-03),em_real_k8(-1.136887572032341D-03) & ! Venus ,em_real_k8(-1.720173260060367D-02),em_real_k8(-3.148295197112937D-03),em_real_k8(8.961124784916037D-07) & ! Earth ,em_real_k8(-1.198027375283752D-02),em_real_k8(-5.092533276958321D-03),em_real_k8(1.875893458494313D-04) & ! Mars ,em_real_k8(3.174273016845083D-03),em_real_k8(7.161028006264430D-03),em_real_k8(-1.007937605937214D-04) & ! Jupiter ,em_real_k8(-4.520953258457782D-04),em_real_k8(-5.588245332485017D-03),em_real_k8(1.148915369554184D-04) & ! Saturn ,em_real_k8(2.706976434750568D-04),em_real_k8(3.738322243478300D-03),em_real_k8(1.038310207754998D-05) & ! Uranus ,em_real_k8(1.746288049471184D-03),em_real_k8(2.613667805270097D-03),em_real_k8(-9.397429335776478D-05)/ ! Naptune ! frame = 0 time = 0.0D0 p_time = 0.0D0 ! 1000 CONTINUE ! IF(MOD(frame,sample_interval)==0)THEN WRITE(*,'(F8.2,18E21.12)')time,((s(axis,body1),axis = 1,2),body1 = 1,9) ENDIF frame = frame+1 ! Note that time will jitter with finite word length. Therefore: time = REAL(frame)/REAL(frames_per_day) dt = time-p_time p_time = time ! DO body1 = 1,9 DO body2 = body1+1,9 DO axis = 1,3 ! Difference in position ds(axis,body1,body2) = s(axis,body1)-s(axis,body2) ! Square of difference in position s2(axis) = ds(axis,body1,body2)*ds(axis,body1,body2) ENDDO ! Square of radial distance r2(body1,body2) = s2(1)+s2(2)+s2(3) ! Radial distance r(body1,body2) = SQRT(r2(body1,body2)) ! Total gravitation field from body2 gf(body1,body2) = g*mass(body2)/r2(body1,body2) ! Total gravitational field from body1 gf(body2,body1) = g*mass(body1)/r2(body1,body2) ! Cosines and vector fields to and from each body DO axis = 1,3 proj(axis,body1,body2) = ds(axis,body1,body2)/r(body1,body2) proj(axis,body2,body1) = -proj(axis,body1,body2) vf(axis,body1,body2) = proj(axis,body1,body2)*gf(body1,body2) vf(axis,body2,body1) = proj(axis,body2,body1)*gf(body2,body1) ENDDO ENDDO ENDDO ! ! Sum accelerations and integrate DO axis = 1,3 ! Total accelerations of each body a(axis,1) = vf(axis,1,2)+vf(axis,1,3)+vf(axis,1,4)+vf(axis,1,5)+vf(axis,1,6)+vf(axis,1,7)+vf(axis,1,8)+vf(axis,1,9) a(axis,2) = vf(axis,2,1)+vf(axis,2,3)+vf(axis,2,4)+vf(axis,2,5)+vf(axis,2,6)+vf(axis,2,7)+vf(axis,2,8)+vf(axis,2,9) a(axis,3) = vf(axis,3,1)+vf(axis,3,2)+vf(axis,3,4)+vf(axis,3,5)+vf(axis,3,6)+vf(axis,3,7)+vf(axis,3,8)+vf(axis,3,9) a(axis,4) = vf(axis,4,1)+vf(axis,4,2)+vf(axis,4,3)+vf(axis,4,5)+vf(axis,4,6)+vf(axis,4,7)+vf(axis,4,8)+vf(axis,4,9) a(axis,5) = vf(axis,5,1)+vf(axis,5,2)+vf(axis,5,3)+vf(axis,5,4)+vf(axis,5,6)+vf(axis,5,7)+vf(axis,5,8)+vf(axis,5,9) a(axis,6) = vf(axis,6,1)+vf(axis,6,2)+vf(axis,6,3)+vf(axis,6,4)+vf(axis,6,5)+vf(axis,6,7)+vf(axis,6,8)+vf(axis,6,9) a(axis,7) = vf(axis,7,1)+vf(axis,7,2)+vf(axis,7,3)+vf(axis,7,4)+vf(axis,7,5)+vf(axis,7,6)+vf(axis,7,8)+vf(axis,7,9) a(axis,8) = vf(axis,8,1)+vf(axis,8,2)+vf(axis,8,3)+vf(axis,8,4)+vf(axis,8,5)+vf(axis,8,6)+vf(axis,8,7)+vf(axis,8,9) a(axis,9) = vf(axis,9,1)+vf(axis,9,2)+vf(axis,9,3)+vf(axis,9,4)+vf(axis,9,5)+vf(axis,9,6)+vf(axis,9,7)+vf(axis,9,8) ! ! Integrate for each body DO body1 = 1,9 IF(frame==1)THEN ! For the first frame we make the previous acceleration equal to the current. This causes an Euler step for ! the velocity integration in this frame. pa(axis,body1) = a(axis,body1) ENDIF ! Keep previous velocities for Tustin's method pv(axis,body1) = v(axis,body1) ! Trapezoidal integration for velocities v(axis,body1) = v(axis,body1)+(3*a(axis,body1)-pa(axis,body1))*dt*0.5D0 ! Tustin's method for positions s(axis,body1) = s(axis,body1)+0.5D0*(v(axis,body1)+pv(axis,body1))*dt ! Keep previous accelerations for next frame pa(axis,body1) = a(axis,body1) ENDDO ENDDO ! end_run = frame>=end_frame IF(.NOT.end_run)GOTO 1000 ! WRITE(*,'(/,F8.2,18E21.12)')time,((s(axis,body1),axis = 1,2),body1 = 1,9) ! ! END PROGRAM