Physical time (years) in simulations

[joseph2015]

Hello

My name is Joseph I am a physics student, I am new to Fortran, I know how to write a code for solving differential equations, but now I have a physics equation for evolving star rotation as a function of mass. and I do not know how to evolve (using any numerical method, such as Euler) so that my code progress and show physical time like mass and rotation in x and y-axis and evolution in years.
any clues, or sample code or useful links please?

Regards.

Joseph

 

[DrClaude]

So long as you are using units consistently, there should be no difficulty. If you are using seconds (for instance, if velocities are in m/s), then you simply need to divide the time by the number of seconds in a year.

 

[jedishrfu]

While this doesn't answer your question directly., I thought it was pretty cool and useful for more serious simulation development in Fortran:

https://github.com/jacobwilliams/Fortran-Astrodynamics-Toolkit

Here's an n-body simulation written in a variety of languages including Fortran that may help you figure out what you need to do:

https://rosettacode.org/wiki/N-body_problem

Lastly, here's a more understandable Java implementation of N-body where these methods actually do the integration of the equations of motion:

        private void computeAccelerations() {
            for (int i = 0; i < bodies; ++i) {
                accelerations[i] = origin;
                for (int j = 0; j < bodies; ++j) {
                    if (i != j) {
                        double temp = gc * masses[j] / Math.pow((positions[i].minus(positions[j])).mod(), 3);
                        accelerations[i] = accelerations[i].plus(positions[j].minus(positions[i]).times(temp));
                    }
                }
            }
        }
        private void computeVelocities() {
            for (int i = 0; i < bodies; ++i) {
                velocities[i] = velocities[i].plus(accelerations[i]);
            }
        }
        private void computePositions() {
            for (int i = 0; i < bodies; ++i) {
                positions[i] = positions[i].plus(velocities[i]).plus(accelerations[i].times(0.5));
            }
        }
        public void simulate() {
            computeAccelerations();
            computePositions();
            computeVelocities();
            resolveCollisions();

https://rosettacode.org/wiki/N-body_problem#Java

Notice they compute the accelerations first then the velocities and finally the positions of the bodies. The simulation repeats this sequence over and over. The computeAccelerations() actually contains the differential equations of motion embedded in the code.

 

[joseph2015]

Thank you very much, that helps alot.

 

[Svein]

The Network Time Protocol (https://en.wikipedia.org/wiki/Network_Time_Protocol) defines a time format and discusses how to transmit time information. The only trouble is that it is defined in "C" terminology. Part of one of my NTP header files follows:

/*---------------------------------------------------*/
/* This is a header file for an SNTP implementation. */
/*                                                   */
/* Start by defining an SNTP timestamp:              */
/*---------------------------------------------------*/

typedef struct _TS  {
  long  TSseconds;
/* -- It should, of course, have been "unsigned long", but -- */
/* -- doing it this way eases the calculations             -- */
  unsigned long  TSfraction;
  } timeStamp;

/*----------------------------------------------------------------------------------*/
/* Below is a description of the NTP/SNTP Version 4 message format, which follows   */
/* the IP and UDP headers. This format is identical to that described in RFC-1305,  */
/* with the exception of the contents of the reference identifier field.            */
/* The header fields are defined as follows:                                        */
/*                                                                                  */
/*                            1                   2                   3             */
/*        0 1 2 3 4 5 6 7 8 9 0 1 2 3 4 5 6 7 8 9 0 1 2 3 4 5 6 7 8 9 0 1           */
/*       +-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+          */
/*       |LI | VN  |Mode |    Stratum    |     Poll      |   Precision   |          */
/*       +-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+          */
/*       |                          Root Delay                           |          */
/*       +-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+          */
/*       |                       Root Dispersion                         |          */
/*       +-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+          */
/*       |                     Reference Identifier                      |          */
/*       +-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+          */
/*       |                                                               |          */
/*       |                   Reference Timestamp (64)                    |          */
/*       |                                                               |          */
/*       +-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+          */
/*       |                                                               |          */
/*       |                   Originate Timestamp (64)                    |          */
/*       |                                                               |          */
/*       +-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+          */
/*       |                                                               |          */
/*       |                    Receive Timestamp (64)                     |          */
/*       |                                                               |          */
/*       +-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+          */
/*       |                                                               |          */
/*       |                    Transmit Timestamp (64)                    |          */
/*       |                                                               |          */
/*       +-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+          */
/*       |                 Key Identifier (optional) (32)                |          */
/*       +-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+          */
/*       |                                                               |          */
/*       |                                                               |          */
/*       |                 Message Digest (optional) (128)               |          */
/*       |                                                               |          */
/*       |                                                               |          */
/*       +-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+          */
/*----------------------------------------------------------------------------------*/

struct NTPmsg  {
  unsigned char li_vn_mode;  /* -- LI | VN | Mode | -- */
  unsigned char Stratum;
  char Poll;
  char Precision;
  long Root_Delay;
  long Root_Dispersion;
  unsigned long refID; /* -- Reference Identifier -- */
  timeStamp refTS;    /* -- Reference Timestamp (64) -- */
  timeStamp origTS;   /* -- Originate Timestamp (64) -- */
  timeStamp rxTS;     /* -- Receive Timestamp (64)   -- */
  timeStamp txTS;     /* -- Transmit Timestamp (64)  -- */
  };

#define SNTP_BUF_SIZE sizeof(struct NTPmsg);

 

[Svein]

Clarifying my last post: The timestamp format is sort of universal, therefore it is usually a good idea to stick with it. There are some caveats, however:

• The specified timestamp rolls over every 232 seconds (136 years) and has a theoretical resolution of 2⁻³² seconds (233 picoseconds). NTP uses an epoch of January 1, 1900 so the first rollover will be on February 7, 2036
• The timestamps used by GPS and IEEE1588 have a similar format, but they differ in details. The rollover problem is the same, however.
• NTPv4 introduces a 128-bit time format: 64 bits for the second and 64 bits for the fractional-second. This format may be of more use to you since "the 64-bit second value is enough to provide unambiguous time representation until the universe goes dim."
• Unfortunately, I do not know how to specify the number of bits in an integer value in FORTRAN.

 

[jedishrfu]

Some info on Fortran integer values:

http://www-classes.usc.edu/engr/ce/108/text/fbk01.htm

http://earth.uni-muenster.de/~joergs/doc/f90/unix-um/dfum_034.html

Basically, the datatype specification does it:

integer*1 for 8 bits (aka byte)
integer*2 for 16 bits (aka short)
integer*4 for 32 bits
integer*8 for 64 bits (aka long)
integer*16 for 128 bits (dont think this exists in today's fortran?)

It looks like you'd have to store the UTC as two integer*8 values.

 

[joseph2015]

Some info on Fortran integer values:

http://www-classes.usc.edu/engr/ce/108/text/fbk01.htm

http://earth.uni-muenster.de/~joergs/doc/f90/unix-um/dfum_034.html

Basically, the datatype specification does it:

integer*1 for 8 bits (aka byte)
integer*2 for 16 bits (aka short)
integer*4 for 32 bits
integer*8 for 64 bits (aka long)
integer*16 for 128 bits (dont think this exists in today's fortran?)

It looks like you'd have to store the UTC as two integer*8 values.

what these have to do with my question? I am really lost!

 

[jedishrfu]

We don’t know what level of detail you need. When folks saw UTC they thought time standards.

What I saw in your post was simply understanding how to setup a numerical computation. I gave you references that should help.

Later, if you need to be precise with the time then perhaps @Svein and my last post will help you extend the time variable to be utc compliant.

 

[joseph2015]

We don’t know what level of detail you need. When folks saw UTC they thought time standards.

What I saw in your post was simply understanding how to setup a numerical computation. I gave you references that should help.

Later, if you need to be precise with the time then perhaps @Svein and my last post will help you extend the time variable to be utc compliant.

thank you very much.