Milankovitch Orbital Data Viewer

Select start and end of timeseries

Compare Distribution of Insolation at Two Different Times

Sensitivity of Insolation to Different Orbital Parameters

Select orbital parameters, compare to present-day

Milankovitch Orbital Parameters over Time

Based on the calculations of:
J. Laskar et al (2004) A long term numerical solution for the insolation quantities of the Earth. Astronomy and Astrophysics 428, 261-285, DOI: 10.1051/0004-6361:20041335
Go to Celestial Mechanics at Observatoire de Paris
Download paper Laskar et al (2011)
download orbital data

Global and Seasonal Distribution of Incoming Solar Radiation

Compare distribution of daily mean radiation at two different times:

Global/Seasonal Distribution of Differences

Difference: Perturbed Orbit - Reference Time:

Gravitational interactions among the many bodies in the Solar System continuously induce small perturbations in the Earth's orbit around the Sun, which in turn cause subtle changes in the distribution of solar radiation by latitude and season over geologic time. In the early 20th Century, Serbian mathematician Milutin Milankovitch proposed that these variations might be related to the timing of the Ice Ages. Because these changes are pseudo-cyclical, they are known as Milankovitch Cycles. For climate, the most import variations are changes in the eccentricity of the ellipse that describes the Earth's orbit, the tilt of the Earth's axis, and the precession of the equinoxes) which affects the seasonal distribution of radiation.

These changes can be predicted (and hindcast) by numerical integration of the multibody equations of celestial mechanics. Jacques Laskar and his colleagues at l'Observatoire de Paris has performed these calculations for 250 million years into the past and for 20 million years into the future. The data and plots here are derived from the output of their calculations, which they have generously made available on their website. The full model, calculations, results, and analysis is described in their paper (Laskar et al, 2011; available from Astronomy and Astrophysics via Open Access ).

The daily mean top-of-the atmosphere incoming solar radiation (insolation) by latitude and season is determined from the orbital parameters using well-known formulae based on spherical trigonometry (see derivation on wikipedia). For all the details, please see the “website code” tab to the right.

Some excellent web-based educational resources on Milankovitch forcing of climate are:

How this website works (including all the code!)

This website is controlled using the R package “shiny.” This one was a lot of work!

There are 10 components:

  • A script (laskar.R) that reads the orbitral data from a file
  • A function (Qday.65N.solstice.R) that calculates daily mean insolation at 65N from the given orbital parameters
  • A script plot.timeseries.R that draws the nice plots on the browser
  • A function global.insol.R that computes daily mean top-of-atmosphere insolation by latitude & season from Earth's orbital parameters (eccentricity, obliquity, and longitude of perihelion)
  • A script that plots the distribution of insolation by latitude & season
  • A script compare.two.times that compares the distributon of insolation at two different times chosen by the user
  • A script sensitivity.Qday that lets the user explore the sensitivity of the distribution of insolation to the Earth's orbital parameters
  • a function that write a text file of orbital data and lets the user download it through the web browser
  • The user interface (ui.R) that controls the sliders and passes inputs to the model
  • The “server” (server.R) that responds to user interface events in the web browser

Scroll down or click links in the list above to read all about it!

Read and subset orbital data between user-specified times


laskar <- function(start.time=-1000, end.time=20){

  # Read timeseries of Earth's orbital eccentracity, obliquity, and precession
  # from a big binary file based on the data of Laskar et al (2004)

  # Subset the data between start.time and end.time (thousands of years from now)
  # then calculate top-of-atmosphere daily mean insolation at 65N June solstice


  data <- subset(orbit, kiloyear>= start.time & kiloyear<= end.time)
  insolation <- Qday.65N.solstice(eccentricity=data$eccentricity, 
  return(data.frame(time=data$kiloyear, eccentricity=data$eccentricity,
              obliquity=data$obliquity, perihelion=data$perihelion,

Calculate top-of-atmosphere insolation at 65N at June Solstice


Qday.65N.solstice <- function(S0=1367, eccentricity=0.0167,
                              obliquity=0.4091, perihelion=1.7963){

  # Given Earth's orbital parameters (eccentricity, obliquity, & perhihelion),
  # Calculate the daily mean top-of-atmosphere insolation at 65N June solstice

  # For a derivation see

  latitude <- 65. * pi / 180. <- S0*(1+eccentricity*sin(perihelion+pi))^2 *sin(latitude)*sin(obliquity) 


Make a 4-panel “ladder plot” of orbital parameters over time


plot.timeseries <- function(orbit){

  # Make a "ladder" plot of orbital parameters over time 

  # attach links names in the list orbit to variables in the script

  # Remember original plotting parameters and then set them the way we want
  orig.par <- par(no.readonly=TRUE)
  par(mfrow=c(4,1), mar=c(0,4,2,1), oma=c(5,1,0,0), cex.lab=1.3, cex.axis=1.3)

  # Top plot is eccentricity of Earth's orbit
       xaxt='n', xlab='',ylab='')
  grid(col='black', ny=0)
  abline(v=0, lwd=2)

  # Second plot is obliquity of rotational axis
       xaxt='n', ylab='Degrees',xlab='')
  grid(col='black', ny=0)
  abline(v=0, lwd=2)

  # Third plot is the precession of the equinox around the orbit
       main=expression(bold(paste('Precession Index:   ',e(1+sin(bar(omega)))))),
                  xaxt='n', xlab='',ylab='')
  grid(col='black', ny=0)
  abline(v=0, lwd=2)

  # Bottom plot is the resulting changes in summer solar radiation at 65N 
       main='Mean Daily Insolation at 65N on Summer Solstice', 
       ylab='Watts/m2', sub='kiloyears')
  grid(col='black', ny=0)
  abline(v=0, lwd=2)

  # Add an overall x-axis to the page
  mtext(side=1, line=2.5, 'kiloyears from present')

  # reset plotting parameters to original values

  # release variable names 


Calculate distribution of insolation by latitude and season


global.insol <- function(S0=1367, ecc=0.0167, eps=0.4091, perih=1.7963,
                         nLats=181, nLons=181){

  # Given Earth's orbital parameters (eccentricity, obliquity, and precession),
  # calculate daily mean top-of-atmosphere incoming solar radiation 
  # for each latitude and time of year, then return the results as a 2D array

  # (longitude here refers to time of year ... 
  #  it's the "longitude" of the Earth in it's orbit around the Sun:
  #     Lon=0 is the March Equinox, Lon=90 is June Solstice
  #     Lon=180 is September Equinox, Lon=270 is December Solstice)

  # Most of this function is just prep, the real calculation happens 
  # just before the return value is assigned: Qday as a function of orbit

  # For details and a derivation see

  # Set up arrays to hold lat, lon, and insolation for each point
  lat <- replicate(nLats,NA)
  lon <- replicate(nLons,NA)
  Qday <- matrix(nrow=nLats, ncol=nLons)

  # Loop over all the points 
  for (j in 1:nLats) for (i in 1:nLons) {

    # latitude and longitude for this [i,j] pair
    lat[j] <- (-90 + (j-1) * 180/(nLats-1)) * pi / 180
    lon[i] <- (-180 + (i-1) * 360/(nLons-1)) * pi / 180

    # Declination for this time of year (longitude)
    # declination is the latitude where the Sun is directly overhead
    dec <- eps * sin(lon[i])

    # Hour angle at sunrise
    arg <- tan(lat[j])*tan(dec)
    if (arg > 1) {
      h0 <- pi # perpetual summer daylight
    } else if (arg < -1) {
      h0 = 0} # perpetual winter darkness
    else {
      h0 <- acos(-arg)  # normal sunrise/sunset

    # Earth-Sun distance (ratio of mean to actual)
    R0.div.RE <- 1 + ecc * cos(lon[i] - perih)

    # Daily mean insolation for this lat and time
    Qday[i,j] <- S0/pi * (R0.div.RE)^2 * 
      (h0 * sin(lat[j]) * sin(dec) + cos(lat[j]) * cos(dec) * sin(h0))
  return(insol=list(lat=lat*180/pi, lon=lon*180/pi, Qday=Qday))

Plot distribution of insolation by latitude and season <- function(insol, plot.title='Reference time', diff=FALSE){

  # Make a color-filled contour plot of global TOA insolation by lat & season
  # optionally use a nice difference color bar

  library(fields)  # Need this to make the image.plot

  # x-coordinate is longitude of Earth in orbit relative to March Equinox
  # This is essentially the date in the year, so label with month names
  x.labels =c('SEP','OCT','NOV','DEC','JAN','FEB','MAR',

  # y-coordinate is latitude on Earth
  y.labels = c('SP','60 S','30 S','EQ','30N','60 N','NP')

  # Optionally center colors on zero for difference plots
  if (diff) {
    z.limits <- c(-max(abs(insol$Qday)),max(abs(insol$Qday)))
  }  else {
    z.limits <- c(min(insol$Qday),max(insol$Qday))

  # Draw an image using colors to represent the insolation
  image.plot(insol$lon,insol$lat,insol$Qday, xlim=c(-180,180), ylim=c(-90,90),
             zlim=z.limits, xlab='', ylab='Latitude', axes=FALSE)

  # Add contours lines for insolation 
  contour(insol$lon,insol$lat,insol$Qday, axes=FALSE, add=TRUE)

  # Add nicely-labeled axes for months and latitude
  axis(1, at=seq(-180,180,30), labels=x.labels)
  axis(2, at=seq(-90,90,30), labels=y.labels)

  # Add titles 
  title(expression(W ~~ m^-2), line=1)


Compare distribution of insolation at two specified times


compare.two.times <- function(orbit, time.1, time.2) {

  # Given a timeseries of the Earth's orbital parameters, select two times
  # (specified in input as time.1 and time.2, thousands of years from now)

  # For each specified time, calculate and plot the distribution of
  # insolation at the top of the atmoshere by latitude and season

  # Make a nice color-filled plot of geography of sunlight at time.1
  # and subtract to show differences in distribution of insolation at time.2

  # Extract orbital parameters for time.1
  first <- subset(orbit, time == time.1)

  # Extract orbital parameters for time.2
  second <- subset(orbit, time == time.2)

  # Calculate geography of top-of-atmosphere sunshine for time.1
  first.insol <- global.insol(ecc=first$eccentricity,

  # Calculate geography of top-of-atmosphere sunshine for time.2
  second.insol <- global.insol(ecc=second$eccentricity,

  # Subtract insolation distribution for time.2 from time.1
  insol.diff <- list(lat=first.insol$lat, 
                     Qday=first.insol$Qday - second.insol$Qday)

  # Set graphics parameters to show a vertical stack of two plots on the page
  orig.par <- par(no.readonly=TRUE)

  # Plot the distribution of insolation at time.1 by latitude & season, plot.title=paste(time.1,'kyr'))

  # Plot the distribution of insolation differences at time.1, diff=TRUE,
                    plot.title=paste(time.1,'kyr minus',time.2,'kyr'))

  # Reset graphics parameters to their original values

Calculate sensitivity of insolation by latitude & season to orbital parameters


sensitivity.Qday <- function(eccentricity, obliquity, perihelion) {

  # Given the Earth's orbital parameters, calculate and plot the distribution of
  # insolation at the top of the atmoshere by latitude and season, and compare 
  # to the distribution in the present day

  # Calculate present-day distribution of insolation by latitude & season <- global.insol()

  # Calculate perturbed distribution of insolation by latitude & season
  perturbation <- global.insol(ecc=eccentricity,

  # Subtract perturbed from present-day insolation by latitude & season 
  insol.diff <- list($lat,$lon,
           $Qday - perturbation$Qday)

  # Set up graphics parameters for a column of two plots stacked on the page
  orig.par <- par(no.readonly=TRUE)

  # Draw the plot of present-day insolation, plot.title='Present Day Conditions')

  # Draw the plot of differences between perturbed and present conditions, diff=TRUE,
                    plot.title='Differences with Perturbed Orbital Parameters')

  # Label the lower plot with the perturbed orbital parameters
  mtext(substitute(e==esub ~","~ epsilon==osub^o ~","~ omega==psub^o,
                        osub=obliquity * 180/pi,
                        psub=perihelion* 180/pi)),
        side=1, line=3, cex=2)

  # Reset graphics parameters to original values

Write orbital data to a text files & download via web browser <- function(file, {
  cat('Data from\n')
  cat('Five columns for each time as follows:\n')
  cat('  time: thousands of years from present\n')
  cat('  eccentricity: (aphelion-perihelion)/(aphelion+perihelion)\n')
  cat('  obliquity: radians\n')
  cat('  perihelion: angle between equinox & perihelion, radians\n')
  cat('  insolation: daily mean TOA solar flux at 65N summer solstice, W/m2\n')
  write.table(, file=file,row.names=F,col.names=T,
              quote=FALSE, append=TRUE, sep="\t")

This is the user interface that actually controls this website



# Define user interface for Milankovitch Data Viewer


  #  Application title
  headerPanel("Milankovitch Orbital Data Viewer"),

  # Sidebar on the left with controls for the user to specify times

    # Display an image that shows the definition of each orbital parameter
    img(src='orbital.jpg', align='left'),

    # Allow the user to choose the start time and duration of the orbital data
    conditionalPanel(condition="input.tabset == 'tab1'", wellPanel(
      h4('Select start and end of timeseries'),            
      sliderInput("start.year", "Beginning of Timeseries (kiloyears from now)", 
                  min = -20000, max = 10000, value = -1000, step= 10),
      sliderInput("duration", "Duration: Thousands of Years", 
                  min = 1, max = 10000, value = 1050, step= 10)

    # Choose two times to compare geography of insolation between them
    conditionalPanel(condition="input.tabset == 'tab2'", wellPanel(
      h4('Compare Distribution of Insolation at Two Different Times'),            
      sliderInput("first.time", "First time (kiloyears from now)", 
                  min = -1000, max = 1000, value = 0, step= 10),
      sliderInput("second.time", "Second time (kyr from now)", 
                  min = -1000, max = 1000, value = -100, step= 10)

    # Choose arbitrary orbital parameters and compare insolation to present
    conditionalPanel(condition="input.tabset == 'tab3'", wellPanel(
      h4('Sensitivity of Insolation to Different Orbital Parameters'),
      h5('Select orbital parameters, compare to present-day'),
      sliderInput('eccentricity', 'Eccentricity',
                  min=0, max=0.060, value=0.0167, step=0.001),
      sliderInput('obliquity', 'Obliquity (degrees)',
                  min=22, max=24.5, value=23.44, step=0.01),
      sliderInput('long.perihelion', 'Longtidue of Perihelion (omega)',
                  min=-180, max=180, value=70.3, step=0.1)

  # Show a tabset in the main panel of the browser that displays model output

      # First tab shows plots of model output for both control and modified atmosphere
      # Also displays explanatory text, links to sources, and offers data download          
      tabPanel("Timeseries", value='tab1',
               h4('Milankovitch Orbital Parameters over Time'), 

               # Explanatory text about the data and sources
               h5("Based on the calculations of:"),
               helpText("J. Laskar et al (2004) A long term numerical solution for the 
               insolation quantities of the Earth. Astronomy and Astrophysics 
               428, 261-285, DOI: 10.1051/0004-6361:20041335"),
               h5('Go to ',a('Celestial Mechanics at Observatoire de Paris',
               h5('Download paper', a('Laskar et al (2011)',

               # Offer to download orbital data through user's web browser
               downloadButton('downloadData',label='download orbital data'),

               # Make the nice 4-panel tiemseries plot
               tags$style(type="text/css", ".tab-content { overflow: visible; }")), 

      # Second tab shows global/seasonal distribution of insolation at reference time
      tabPanel("Compare", value='tab2',
               h4('Global and Seasonal Distribution of Incoming Solar Radiation'), 
               h5('Compare distribution of daily mean radiation at two different times:'),
               tags$style(type="text/css", ".tab-content { overflow: visible; }")), 

      # Third tab compares insolation distribution between ref time and perturbations
      tabPanel("Sensitivity", value='tab3',
               h4('Global/Seasonal Distribution of Differences'), 
               h5('Difference: Perturbed Orbit - Reference Time:'),
               tags$style(type="text/css", ".tab-content { overflow: visible; }")), 

      # Fourth tab displays a brief model description
      tabPanel("Description", value='tab4',

      # Fifth tab displays website code
      tabPanel("Website code", value='tab5',

This is the server function that responds to user input



# Source required R scripts

shinyServer(function(input, output) {

# Capture orbital timeseries for later display <- reactive({

  # 4-panel plot of orbital timeseries including insolation at 65N
  output$timeseries <- renderPlot(
    plot.timeseries(, height=800, width=600)

  # Download orbital data as a text file from web browser
  output$downloadData <- downloadHandler(
    filename = '',
    content = function(file) {,

  # Contour plot of insolation differences for given orbital parameters
  output$compare.times <- renderPlot(
      compare.two.times(,input$first.time, input$second.time), 

  # Explore sensitivity to orbital parameters:
  output$sensitivity <- renderPlot(
    sensitivity.Qday(input$eccentricity, input$obliquity * pi/180,
                     input$long.perihelion * pi/180),
    height=950, width=600)