%GCL
// -----------------------------------------------------------------------------
// Plot La2004 data
// la2004_plot100.gcl - created 2006-Jan-23, last updated 2006-Jan-29
// (c) Erik De Sonville, http://users.skynet.be/erik.desonville/
// Scientific advisor (geology): Dr. Thierry Moorkens
// Script language: disgcl (H. Michels), http://www.dislin.de/
// -----------------------------------------------------------------------------
// Input file: insoln.la2004.btl.100.asc
//   This is an orbital-parameters file from the La2004 model,
//   namely the file with single-precision data covering 0 to -101 Myr
// Output file:
//   .png plot of orbital parameter (as defined/selected)
// -----------------------------------------------------------------------------
// Script status: experimental (quick and not very clean)
// - parameterization of graphs
// - various plots via parameter patches
// - later: interactive widgets for user selections
// -----------------------------------------------------------------------------

//*** parameter patch zone *****************************************************
  swv_of = 0
  swv_on = 1
//--> define toggle switches here
sw_debugm = swv_of	// show debug messages
sw_debugn = swv_of	// show debug numerics
sw_plot   = swv_on	// generate plot

//--> define graph subtitle here
tit2 = 'Demo graph, prepared 2006-Jan-29'

  swv_ecc=1	// eccentricity
  swv_obl=2	// obliquity
  swv_pri=3	// precession index
//--> select orbital parameter here
sel_btl=swv_pri
  if      (sel_btl==swv_ecc)
    btlstr='eccy'
  else if (sel_btl==swv_obl)
    btlstr='obly'
  else if (sel_btl==swv_pri)
    btlstr='priy'
  else
   goto einp
  end if
  
  swv_mio=1	// Miocene, Holbourn
  swv_ru1=2	// Boom clay
  swv_ru2=3 // Boom clay ext.
  swv_yp1=4	// Egem sand
  swv_yp2=5 // Roubaix clay
  swv_yp3=6	// Knokke clay
//--> select geological stage here
sel_age=swv_yp2
  if      (sel_age==swv_mio)
    perstr = '-13.7-14.2'
    yr0 = -13700
    per =    500
  else if (sel_age==swv_ru1)
    perstr = '-28.0-31.0'
    yr0 = -28000
    per =   3000
  else if (sel_age==swv_ru2)
    perstr = '-31.0-34.0'
    yr0 = -31000
    per =   3000
  else if (sel_age==swv_yp1)
    perstr = '-50.5-51.5'
    yr0 = -50500
    per =   1000
  else if (sel_age==swv_yp2)
    perstr = '-52.0-53.5'
    yr0 = -52000
    per =   1500
  else if (sel_age==swv_yp3)
    perstr = '-54.5-55.5'
    yr0 = -54500
    per =   1000
  else
    goto einp
  end if

  outfile = '.\data\astro_rupel_plots\' + btlstr + perstr + '.png'
  if (sw_debugm==swv_on) print 'outfile:',outfile
//*** parameter patch zone end *************************************************

// Variables

pi = 3.14159265359
todeg = 180./pi

stat = 0            // return status of functions

yr9 = yr0-per
if (yr9>=0) goto einp
if (yr0<-101000) goto einp
if (yr9>=yr0) goto einp

nsk  = -yr0
nyr  = -yr9+yr0+1

eccy = falloc(nyr)  // orbital element arrays
obly = falloc(nyr)
priy = falloc(nyr)

hbtl = -1		// file handle

// Read data into arrays

hbtl = fopen('.\data\btl\insoln.la2004.btl.100.asc','r')
  if (hbtl==-1) goto ebtl

  // skip the non-relevant part
do i=0,nsk-1
  stat = fscanf(hbtl,'%f%f%f%f',kyear,eccen,obliq,lonph)
    if (stat==-1) goto ebtl
end do
if (sw_debugm==swv_on) print 'skip done'

  // read arrays

do i=0,nyr-1
  stat = fscanf(hbtl,'%f%f%f%f',kyear,eccen,obliq,lonph)
    if (stat==-1) goto ebtl
  eccy[i] = eccen
  obly[i] = obliq
  priy[i] = eccen * sin(lonph)
end do
stat = fclose(hbtl)
  if (stat==-1) goto ebtl
if (sw_debugm==swv_on) print 'btl read'
if (sw_debugn==swv_on) print 'first/last eccy',eccy[0],eccy[nyr-1]
if (sw_debugn==swv_on) print 'first/last obly',obly[0],obly[nyr-1]
if (sw_debugn==swv_on) print 'first/last priy',priy[0],priy[nyr-1]

// Prepare arrays for plotting

x0 = float(-yr0)
x9 = float(-yr9)

if      (per==500)
  xs = 100.
  xt = 5
else if (per==1000)
  xs = 200.
  xt = 4
else if (per==1500)
  xs = 500.
  xt = 5
else if (per==2000)
  xs = 500.
  xt = 5
else if (per==2500)
  xs = 500.
  xt = 5
else if (per==3000)
  xs = 500.
  xt = 5
else
  xs = per/5
  xt = 5
end if

xarr = falloc(nyr)
xarr = xarr+x0
xlabdig=-1
x0 = x0/1000.
x9 = x9/1000.
xs = xs/1000.
xarr = xarr/1000.
xlabdig=1

if (sel_btl==swv_ecc)
  tit1 = 'Eccentricity (from La2004 model)'
  y0 = 0.
  y9 = 0.07
  ys = 0.01
  yname = 'Eccentricity'
  ylabdig = 2
else if (sel_btl==swv_obl)
  tit1 = 'Obliquity (from La2004 model)'
  y0 = 0.38
  y9 = 0.43
  ys = 0.01
  yname = 'Obliquity (rad)'
  ylabdig= 2
else if (sel_btl==swv_pri)
  tit1 = 'Precession Index (calculated from La2004 model)'
  y0 = -0.08
  y9 = +0.08
  ys =  0.02
  yname = 'Precession Index'
  ylabdig = 2
else
  goto einp
end if

// Plot

if (sw_plot==swv_on)

metafl ('cons')
disini ()
hwfont ()
pagera ()

axspos (300,320+1200)
axslen (2000,1200)

name   ('Time before J2000.0 (Ma)','x') 
name   (yname, 'y') 
titlin (tit1,1)
titlin (tit2,2)
ticks  (xt, 'x')
labdig (xlabdig,'x')
labdig (ylabdig,'y')

graf   (x0,x9,x0,xs,y0,y9,y0,ys)
title  ()
color  ('white')

if (sel_btl==swv_ecc) curve (xarr,eccy,nyr)
if (sel_btl==swv_obl) curve (xarr,obly,nyr)
if (sel_btl==swv_pri) curve (xarr,priy,nyr)

imgclp (0,0,820,600)
filmod ('delete')
rpng   (outfile)

disfin ()

end if

// Exits

exitok:
  if (sw_debugm==swv_on) print 'Done okay'
  goto endpgm
einp:
  print 'Error inp'
  goto endpgm
ebtl:
  print 'Error btl'
  goto endpgm
endpgm: