# Insolation.ndf 17Feb07 www.BillHowell.ca Warning: Use Courier constant-width font and no-line-wrap to view this file! Otherwise the code doesn't line up and it's hard to read. # Original code from: J. Laskar, P. Robutel, F. Joutel, M. Gastineau, A.C.M. Correia, B. Levrard 2004 "A long term numerical solution for the insolation quantities of the Earth" Astronomy & Astrophysics manuscript no. La 2004 May 23, 2004 http://www.imcce.fr/fr/presentation/equipes/ASD/preprints/prep.2004/La_2004_prep.pdf # Laskar, J., Joutel, F., Boudin, F. 1993 "Orbital, precessional and insolation quantities for the Earth from -20 Myr to + 10Myr" Astron. Astrophys. 270, 522 http://www.imcce.fr/Equipes/ASD/insola/earth/La93/README.TXT #**************************************** # Laskar etal - insola.f ported to the Q'Nial programming language by www.BillHowell.ca * Howell: View this file with constant-width font (eg Courier) and word wrap turned off * Much of the original content and references is retained. * * * Computation of the various insolation functions * * this program uses * insolsub.f : subroutines for the computation * of insolation * prepinsol.f : subroutines for the binary temporary file * prepsub.f : subroutines for computing elements at * specified time * insola.par : parameters for insola.f * nomclimaneg , * nomclimapos : ASCII files supplied * which contains the necessary * orbital and precession quantities * * * * the following parameters are taken from the NAMELIST insola.par * nomclimapos : ASCII file for the elements t,e,eps, psibar on * positive time. * nomclimaneg : ASCII file for the elements t,e,eps, psibar on * negative time. * nominsolbin : temporary binary file for the elements t,e,eps, * cos(psibar), sin(psibar) * (variables interne = nomfich) * pas : stepsize (years) * datefin : final time (million of years) * datedebut : starting time (million of years) * so : solar constant W*m^-2 * * other parameters: * * nbf : 4 (number of variables) * nbel1 : 4 + 1 * nacd : numbers of n-uplets in a record * * (c) Astronomie et Systemes Dynamiques/ institut de Mecanique Celeste (2003) * J. Laskar, F. Joutel, M. Gastineau * * version 0.82 (7/7/1993) * version 0.83 (27/1/93) * version la2001b4 (26/06/2001) Mickael : changement du fichier de parametres * version la2001b9 (12/02/2002) Mickael : ajout dans le fichier de parametres * version la2003 (04/03/2003) Mickael : fichier de parametres * version la2003 (04/09/2003) Mickael : creation du fichier binaires * (05/11/2003) Mickael : changement du common date en ddate * version la2004 (26/05/2004) Mickael : correction qromb (insolsub.f) #**************************************** #********************************************* * Laskar etal - insolsub.ndf * * subroutines for the computation of insolation * J. Laskar, F. Joutel, F. Boudin * (c) Astronomie et Systemes Dynamiques / Institut de Mecanique Celeste (2003) * * 05 / 11 / 2001 Mickael : modification de l'appel de polint dans qromb * (mauvais type fourni) * 22 / 09 / 2003 Mickael : utilisation de e et pibarh au lieu de ak,ah, psi * 05 / 11 / 2003 Mickael : changement du common date en ddate d * 26 / 05 / 2004 Mickael : correction qromb suivant numerical recipes 2nd ed. * * traduction en Q'Nial de www.qnial.com - Michael, Jenkins & team at Queen's University, Canada * par www.BillHowell.ca février 2007 * #********************************************* #*************************************** * Description of routines, Howell Feb07 * * e - excentricité * eps - l'obliquité * hl - longitude vraie , moyvrai(hlm, e, pibar) * hlm - longitude moyenne, vraimoy(hl , e, pibar) * phi - latitude du point sur la terre * pibar - longitude du perihelie !&!comptee a partir de l'equinoxe de la date + pi (repere geocentrique)!&! * wd - longitude vraie du soleil (en radian) compte a partir de l'equinoxe de la date * * * * vraimoy (hl,e,pibar,hlm) - calcule la longitude moyenne (hlm) * moyvrai (hlm,e,pibar,hl) - calcule la longitude vraie (hl) * cwj (wd,e,pibar,eps,phi) - INSOLATION JOURNALIERE D'UN POINT DE LATITUDE DONNEE * wmcal (mois,e,eps,pibarh,phi) - INSOLATION MENSUELLE D'UN POINT DE LATITUDE DONNEE * wjour (date,e,eps,pibarh,phi) - INSOLATION JOURNALIERE D'UN POINT DE LATITUDE DONNEE * wjcal (datecal,e,eps,pibarh,phi) - INSOLATION JOURNALIERE D'UN POINT DE LATITUDE DONNEE * Telinsol (tps) - tableau des éléments orbitaux (k,h,eps,phi) pour la date tps * (Howell: telinsol not called in Laskar program !!&!!!&!!) NOT called by prepinsol.f nor prepsub.f !!!!! * Global_annual_avg_insolation (e) - INSOLATION ANNUELLE MOYENNEE * * * From standard libraries at !&!l'Observatoire de Paris!&! * (c) Astronomie et Systemes Dynamiques/ institut de Mecanique Celeste (2003) * * QROMB (FUNC,A,B,SS) - ROMBERG INTEGRATION * POLINT (XA, YA, N, X) - used by qromb * TRAPZD (FUNC A B N) - used by qromb * * ******************************************* #********************************* # globals defined here, or used in this file only ; e_du := 0.0 ; e_du_v := f3 ; eps_du := 0.0 ; eps_du_v := f3 ; nbfmax := 20 ; nbfmax_v := f3 ; q_dy := 0.0 ; q_dy_v := f3 ; % for rhomberg routines ; q_nmax := 10 ; q_nmax_v := f3 ; % for rhomberg routines ; phi_du := 0.0 ; phi_du_v := f3 ; pibar_du := 0.0 ; pibar_du_v := f3 ; % for rhomberg integration routines ; q_c := q_nmax reshape 0.0 ; q_c_v := f3 ; % for rhomberg routines ; q_d := q_nmax reshape 0.0 ; q_d_v := f3 ; % for rhomberg routines ; q_xa := q_nmax reshape 0.0 ; q_xa_v := f3 ; % for rhomberg routines ; q_ya := nbfmax q_nmax reshape 0.0 ; q_ya_v := f3 ; % for rhomberg routines ; cwj_o := f3 ; civilization_insolation_o := f3 ; Global_annual_avg_insolation_o := f3 ; Isolation_mensuel_o := f3 ; moyvrai_o := f3 ; QROMB_o := f3 ; POLINT_o := f3 ; Telinsol_o := f3 ; TRAPZD_o := f3 ; vraimoy_o := f3 ; wmcal_o := f3 ; wjour_o := f3 ; wjcal_o := f3 ; # EXTERNAL definitions to allow this file to load alone civilisation_period IS external variable ; civilisation_period_x := f3 ; fichout IS external variable ; fichout_x := f3 ; intrp_tseries IS external variable ; intrp_tseries_x := f3 ; Milankovic IS external variable ; Milankovic_x := f3 ; mois IS external variable ; mois_x := f3 ; n_intrpsteps IS external variable ; n_intrpsteps_x := f3 ; phi IS external variable ; phi_x := f3 ; so IS external variable ; so_x := f3 ; interpolation_for_interval IS external operation ; interpolation_for_interval_x := f3 ; # ********************************************* c subroutine vraimoy(hl,e,pibar,hlm) c c calcule la longitude moyenne (hlm) en fonction de la longitude c vraie (hl), de l'excentricite (e) et de la longitude du perihelie c (pibar) c (c) ASD / BDL (1 / 10 / 92) c correction 27 / 1 / 93 # ********************************************* vraimoy IS OP hl e pibar { beta := (1 - (e * e)) power 0.5 ; eM := hl - pibar ; hl - (2 * (e / 2 + (e power 3 / 8) ) * (1 + beta) * sin eM ) - ( (e * e / 4 ) * (0.5 + beta) * sin (2 * eM) ) + ( (e power 3 ) * (1 / 3 + beta) * sin (3 * eM) ) } # ********************************************* c subroutine moyvrai(hlm,e,pibar,hl) c c calcule la longitude vraie (hl) en fonction de la longitude c moyenne (hlm), de l'excentricite (e) et de la longitude du c perihelie (pibar) c (c) ASD / BDL (1 / 10 / 92) # ********************************************* moyvrai IS OP hlm e pibar { eM := hlm - pibar ; hlm + ( ( 2 * e - (e power 3 / 4)) * sin eM ) + ( 5 / 4 * (e * e) * sin (2 * eM ) ) + ( 13 / 12 * (e power 3) * sin (3 * eM ) ) } # ********************************************* c SUBROUTINE cwj(wd,e,pibar,eps,phi,w) c c INSOLATION JOURNALIERE D'UN POINT DE LATITUDE DONNEE c c A L'ENTREE: c wd : la longitude vraie du soleil (en radian) compte a partir c de l'equinoxe de la date c e : excentricite c pibar: longitude du perihelie comptee a partir de l'equinoxe de la c date + pi (repere geocentrique) c eps: l'obliquite c phi: latitude du point sur la terre c c EN SORTIE: c w:l'insolation journaliere d'un point de latitude donnee c (c) ASD / BDL (1 / 10 / 92) # ********************************************* cwj IS OP wd e pibar eps phi { v := wd - pibar ; rho := (1 - (e * e)) / (1 + (e * cos v) ) ; % Distance terre - soleil (rho): ; delta := arcsin (sin eps * sin wd) ; % Declinaison du soleil (delta): ; aux := pi / 2 - abs delta ; IF (( - aux < phi) and (phi < aux)) THEN % LATITUDES DE LEVER ET COUCHER DE SOLEIL: ; cho := - (tan phi) * (tan delta) ; ho := arccos cho ; % Angle horaire du lever et coucher du soleil (ho): ; ( ( ho * (sin phi) * (sin delta) ) + ( (cos phi) * (cos delta) * (sin ho) ) ) * so / (pi * (rho * rho) ) ELSE a1 := pi / 2 - delta ; a2 := pi / 2 + delta ; IF ((phi >= a1) or (phi <= ( - a2)) ) THEN % LATITUDES SANS COUCHER: ; so * (sin phi) * (sin delta) / (rho * rho) ELSE % LATITUDES SANS LEVER: ; IF ((phi <= ( - a1)) OR (phi >= a2) ) THEN 0 ENDIF ENDIF ENDIF } # ********************************************* * Isolation_mensuel # ********************************************* * Isolation_mensuel IS OP hlm { cwj (moyvrai hlm e_du pibar_du) e_du pibar_du eps_du phi_du } # ********************************************* * c SUBROUTINE wjour(date,e,eps,pibarh,phi,w) c c INSOLATION JOURNALIERE D'UN POINT DE LATITUDE DONNEE c c A L'ENTREE: c date : la longitude vraie du soleil (en degres) comptee a partir c de l'equinoxe vrai c e,eps,pibarh: les elements orbitaux pour la date tps c e : excentricite c eps:l'obliquite c pibarh: longitude du perihelie comptee a partir de l'equinoxe de la c date c phi: latitude du point sur la terre c c c ---------------------------------------------------------------------- - c CALCUL D'INSOLATION c ---------------------------------------------------------------------- - c c longitude du perihelie par rapport a l'equinoxe de reference(pitild), c par rapport a l'equinoxe de la date(pibar), excentricite(e): c c c Longitude vraie comptee a partir de l'equinoxe vrai et c anomalie vraie (wd,v): c c c EN SORTIE: c w:l'insolation journaliere d'un point de latitude donnee c (c) ASD / IMC (22 / 09 / 2003) # ********************************************* * wjour IS OP date e eps pibarh phi { pibar := pibarh + pi ; hl := date * pi / 180 ; cwj hl e pibar eps phi } # ********************************************* * c SUBROUTINE wjcal(datecal,e,eps,pibarh,phi,w) c c INSOLATION JOURNALIERE D'UN POINT DE LATITUDE DONNEE c c A L'ENTREE: c datecal: la longitude moyenne (en degres) du soleil comptee a partir c de l'equinoxe de la date c e,eps,pibarh: les elements orbitaux pour la date tps c e : excentricite c eps:l'obliquite c pibarh: longitude du perihelie comptee a partir de l'equinoxe de la c date c phi: latitude du point sur la terre c c c ---------------------------------------------------------------------- - c CALCUL D'INSOLATION c ---------------------------------------------------------------------- - c c longitude du perihelie par rapport a l'equinoxe de reference(pitild), c par rapport a l'equinoxe de la date(pibar), excentricite(e): c c c EN SORTIE: c w:l'insolation journaliere d'un point de latitude donnee c (c) ASD / BDL (1 / 10 / 92) # ********************************************* * wjcal IS OP datecal e eps pibarh phi { pibar := pibarh + pi ; hlm := datecal * pi / 180 ; % calcul de la longitude vraie ; hlm0 := vraimoy 0. e pibar ; % longitude moyenne au 21 mars ; hlm := hlm0 + hlm ; % longitude moyenne a la date ; wd := moyvrai hlm e pibar ; % longitude vraie a la date ; cwj wd e pibar eps phi } # ********************************************* c SUBROUTINE Global_annual_avg_insolation(e,w) c c INSOLATION ANNUELLE MOYENNEE c c A L"ENTREE: c e:excentricite pour une date donnee c c ---------------------------------------------------------------------- - c CALCUL D'INSOLATION (w) c ---------------------------------------------------------------------- - c Carre de la distance terre - soleil moyennee sur une revolution (rho): c c c EN SORTIE: c w:insolation annuelle moyenne c (c) ASD / BDL (1 / 10 / 92) # ********************************************* Global_annual_avg_insolation IS OP e { so / (4 * sqrt(1 - e power 2)) } # ************************************************* c SUBROUTINE QROMB(FUNC,A,B,SS) c c NUMERICAL RECIPES c ROMBERG INTEGRATION # ************************************************* c 05 / 11 / 2001 Mickael : ajout de 0 (au lieu de 0.) pour l'appel de polint c 26 / 05 / 2004 Mickael : correction du test IF (ABS(DSS).LT. en .LE. c correction effectuee dans "numerical recipes, second edition" POLINT IS OP N X { NONLOCAL q_c q_d q_dy q_xa q_y q_ya q_nmax ; % ; NS := 1 ; DIF := ABS(X - q_xa@1) ; FOR I WITH count N DO DIFT := ABS (X - q_xa@I) ; IF DIFT < DIF THEN NS := I ; DIF := DIFT ; ENDIF ; q_C@I := q_ya@I ; q_D@I := q_ya@I ; ENDFOR ; q_y := q_ya@NS ; NS := NS - 1 ; FOR M WITH count (N - 1) DO FOR I WITH count (N - M) DO HO := q_xa@I - X ; HP := q_xa@[I + M] - X ; W := q_c@[I + 1] - q_d@I ; DEN := HO - HP ; IF DEN = 0.0 THEN writescreen ' pause polint - den=0 error where interpolated point is at known point' ; ENDIF ; % what does this do!&! ; DEN := W / DEN ; q_d@I := HP * DEN ; q_c@I := HO * DEN ; ENDFOR ; IF ( (2 * NS) < (N - M) ) THEN q_dy := q_c@[NS + 1] ; ELSE q_dy := q_d@NS ; NS := NS - 1 ; ENDIF ; q_y:= q_y+ q_dy ; ENDFOR ; % return 2 results ; q_y q_dy } TRAPZD IS OP FUNC A B N { % SAVE IT ; % what does IT do!&! intermediate save in file!&! ; IF N = 1 THEN S := 0.5 * (B - A) * ((apply FUNC A) + (apply FUNC B)) ; IT := 1 ; ELSE TNM := IT ; DEL := (B - A) / TNM ; X := A + (0.5 * DEL) ; SUMT := 0. ; FOR J WITH count IT DO SUMT := SUMT + (FUNC X) ; X := X + DEL ; ENDFOR ; S := 0.5 * (S + ( (B - A) * SUMT / TNM) ); IT := 2 * IT ; ENDIF ; S } QROMB IS OP FUNC A B { error_qromb := o ; EPS := 1.e-6 ; JMAX := 20 ; JMAXP := JMAX + 1 ; K := 5 ; KM := 4 ; S := JMAXP reshape 0.0 ; H := S ; H@1 := 1. ; FOR J WITH count JMAX DO S@J := TRAPZD FUNC A B J ; IF J >= K THEN Lx := J - KM ; SS DSS := POLINT H@Lx S@Lx K 0.0 ; IF (ABS DSS <= (EPS * ABS SS) ) THEN error_qromb := l ; ENDIF ; ENDIF ; S@(J + 1) := S@J ; H@(J + 1) := 0.25 * H@J ; ENDFOR ; IF error_qromb = l THEN Name := Readscreen 'QROMB: Too many steps.' ELSE SS ENDIF } # ********************************************* * c SUBROUTINE wmcal(mois,e,eps,pibarh,phi,w) c c INSOLATION MENSUELLE D'UN POINT DE LATITUDE DONNEE c c A L'ENTREE: c mois : le numero du mois (l'annee est divisee en 12 mois de 30 c degres et le numero 3 correspond a fin fevrier - fin mars) c e,eps,pibarh: les elements orbitaux pour la date tps c e : excentricite c eps:l'obliquite c pibarh: longitude du perihelie comptee a partir de l'equinoxe de la c date c phi: latitude du point sur la terre c c c ---------------------------------------------------------------------- - c CALCUL D'INSOLATION c ---------------------------------------------------------------------- - c c longitude du perihelie par rapport a l'equinoxe de reference(pitild), c par rapport a l'equinoxe de la date(pitbar), excentricite(e): c c Longitude moyenne au 21 mars c c c EN SORTIE: c w:l'insolation mensuelle d'un point de latitude donnee c (c) ASD / IMC (22 / 09 / 2003) # ********************************************* * # Howell - this probably isn't useful (including QROMB) forthe very short timeframes that I am looking at wmcal IS OP mois e eps pibarh phi { NONLOCAL phi_du eps_du e_du pibar_du ; % keep these intermediates, as I don't know if conflicts will arise by making e eps puibarh etc global ; phi_du := phi ; eps_du := eps ; e_du := e ; % globals so that QROMB & Isolation_mensuel can work ; pibar_du := pibarh + pi ; % ; hlm0 := vraimoy 0 e pibar_du ; % longitude moyenne (hlm0) au debut de l'annee!&! ; hlm1 := hlm0 + (mois - 4) * pi * 30 / 180 ; % longitude moyenne (hlm1) au debut du mois ; hlm2 := hlm1 + pi * 30 / 180 ; % longitude moyenne (hlm2) au fin du mois ; % calcul de l'insolation moyenne a l'aide de la methode de Romdberg ; (qromb "Isolation_mensuel hlm1 hlm2) / 30 / pi * 180 } civilization_insolation IS OP fichout latitude mois { LOCAL interval_inputs interval_outpts ; % ; % Solanki ; yr := 0 ; spot := 1 ; irr := 2 ; STf := 3 ; % Milankovic ; yr := 0 ; ecc := 1 ; obl := 2 ; cosp := 3 ; sinp := 4 ; prec := 5 ; % interval_parms ; yr := 0 ; ecc := 1 ; obl := 2 ; cosp := 3 ; sinp := 4 ; prec := 5 ; % interval_inputs ; mois := 0 ; ecc := 1 ; obl := 2 ; pr := 3 ; lat := 4 ; % interval_outpts ; yr := 0 ; ecc := 1 ; obl := 2 ; cosp := 3 ; sinp := 4 ; prec := 5 ; irr := 6 ; inraw := 7 ; insol := 8 ; % ; interval_inputs := n_intrpsteps 5 reshape 0.0 ; % initial assignments same for full problem set ; interval_inputs|[,mois] := n_intrpsteps reshape mois ; % mean monthly, eg 7 => 21 june - 20 july ; interval_inputs|[,lat ] := n_intrpsteps reshape lat ; % latitude du point sur la terre ; interval_outpts := n_intrpsteps 9 reshape 0.0 ; % ; fichout_number := open fichout "w ; fichout_number EACHRIGHT writefile ' ' fichout (link 'Averages for the month of: ' string mois) '' 'Year eccentricity obliquity precession cos_p sin_p irradiance insol_raw insolation' ; % ; FOR interval_start_year WITH civilisation_period DO interval_years := interval_start_year + intrp_tseries ; interval_outpts|[, yr ecc obl cosp sinp] := interpolation_for_interval interval_start_year ; % ; ST_start_index := find interval_start_year Solanki|[,yr] ; interval_STfs := Solanki|[ST_start_index + intrp_tseries, STf] ; % ; interval_inputs|[,ecc obl] := interval_outpts|[,ecc obl] ; interval_inputs|[,pr] := EACH arctan (OUTER / interval_outpts|[,sinp] interval_outpts|[,cosp]) ; interval_outpts|[,inraw] := EACH wmcal rows interval_inputs ; interval_outpts|[,insol] := OUTER * interval_outpts|[,inraw] interval_STfs ; % ; EACH writefile EACH string rows interval_outpts ; ENDFOR ; close fichout_number ; } # enddoc # Main.ndf 17Feb07 www.BillHowell.ca Warning: Use Courier constant-width font and no-line-wrap to view this file! Otherwise the code doesn't line up and it's hard to read. # Original code from: J. Laskar, P. Robutel, F. Joutel, M. Gastineau, A.C.M. Correia, B. Levrard 2004 "A long term numerical solution for the insolation quantities of the Earth" Astronomy & Astrophysics manuscript no. La 2004 May 23, 2004 http://www.imcce.fr/fr/presentation/equipes/ASD/preprints/prep.2004/La_2004_prep.pdf # Laskar, J., Joutel, F., Boudin, F. 1993 "Orbital, precessional and insolation quantities for the Earth from -20 Myr to + 10Myr" Astron. Astrophys. 270, 522 http://www.imcce.fr/Equipes/ASD/insola/earth/La93/README.TXT #**************************************** * * traduction en Q'Nial de www.qnial.com - * Michael, Jenkins & team at Queen's University, Canada * par www.BillHowell.ca février 2007 * * WARNING: This a draft port from fortran to Q'Nial. Errors are guaranteed * * USER INPUTS: Rather than use a binary file for database access to intermediate results, * I've simply used an array in memory (my current focus is limited to the Holocene period). * I'm not sure if an array of 100 k records of 5 to 10 reals each would be too * much memory for this version of Q'Nial, but I think it would be OK. * * Load definitions - with Q'Nial, it is the problem spec file that will drive the actions * (Howell - I prefer not to waste my time with user input prompts, although that option * will eventually be supported) * * Insolation operators (cwj, wmcal, wjour, wjcal, Global_annual_avg_insolation) return a real, * they don't assign to an argument as in theoriginal Laskar etal * * Missing (some other time): * help files with search - in pdf format and pretty fonts etc etc (yuuchh, but that's what people want) * wide range of examples for different types of research * * * Other comments for Q'Nial speedup * replace "x power 2" with "x * x" * don't keep recreating arrays !!! -> QROMB POLINT * * * Bill Howell * #********************************* # globals defined here, or used in this file only ; datedebut := 0.0 ; datedebut_v := f1 ; datefin := 0.0 ; datefin_v := f1 ; itdebut := 0.0 ; itdebut_v := f1 ; % Nombre de lignes lues ; aux := 0.0 ; aux_v := f1 ; npt := 0.0 ; npt_v := f1 ; % Milankovic_rows !! ; civilisation_period := 0.0 ; civilisation_period_v := f1 ; Milankovic_delta_t := 0.0 ; Milankovic_delta_t_v := f1 ; Solanki_delta_t := 0.0 ; Solanki_delta_t_v := f1 ; Complete_inputs_e := f1 ; # EXTERNAL definitions to allow this file to load alone date_conversion IS external variable ; date_conversion_x := f1 ; debut IS external variable ; debut_x := f1 ; interpolation_setup IS external expression ; interpolation_setup_x := f1 ; Milankovic IS external variable ; Milankovic_x := f1 ; Milankovic_rows IS external variable ; Milankovic_rows_x := f1 ; n_intrpsteps IS external variable ; n_intrpsteps_x := f1 ; #**************************************** # Set up files, vérification du pas Complete_inputs IS { NONLOCAL civilisation_period datedebut datefin itdebut npt ; % ; % Milankovic columns ; yr_col := 0 ; ecc_col := 1 ; obl_col := 2 ; cosp_col := 3 ; sinp_col := 4 ; prec_col := 5 ; input_error := l ; % Extra care to make user-define variables are all there!!! ; symbol_list := EACH first symbols 0 ; IF not in "MILANKOVIC symbol_list THEN writescreen '?Error in file: Civilization Milankovic precess, obliq, eccent.txt: ' ; writescreen 'Milankovic data hasnt been loaded' ; input_error := l ; ENDIF ; IF not in "DATE_CONVERSION symbol_list THEN writescreen '?Error in file: Civilization Milankovic precess, obliq, eccent.txt: ' ; writescreen 'Missing date_coversion to change first column of data input to years from ky or My' ; input_error := l ; ENDIF ; IF not in "PAS symbol_list THEN writescreen '?Error in file: Civilization Milankovic precess, obliq, eccent.txt: ' ; writescreen 'Missing stepsize (years - NOT millions of years)' ; input_error := l ; ENDIF ; % data completion & testing continued ; Milankovic|[, yr_col ] := Milankovic|[ , yr_col ] * date_conversion + 2000 ; Milankovic|[, cosp_col] := cos Milankovic|[ , prec_col] ; Milankovic|[, sinp_col] := sin Milankovic|[ , prec_col] ; datedebut := Milankovic@[0, yr_col ] ; % starting time (years) ; datefin := Milankovic@[Milankovic_rows - 1,0] ; % final time (years) ; % Fancy data checks ; % There must be constant time intervals between data dates!! ; IF ~= (Milankovic_delta_t := - [rest, front] Milankovic|[,0]) THEN writescreen '?Error in file: Milankovic - civilization.ndf' ; writescreen 'Milankovic data is provided for an irregular date series - there must be constant years difference between data in the series' ; input_error := l ; ELSE Milankovic_delta_t := first Milankovic_delta_t ; ENDIF ; IF Milankovic_delta_t < 1 THEN writescreen '?Error in file: Milankovic - civilization.ndf ' ; writescreen 'Step size (Milankovic_delta_t) must be >= 1 year' ; input_error := l ; ENDIF ; IF NOT < (datedebut datefin) THEN writescreen '?Error in file: Civilization Milankovic precess, obliq, eccent.txt: ' ; writescreen 'Milankovic dates must start with earliest and end with last' ; input_error := l ; ENDIF ; IF ~= (Solanki_delta_t := - [rest, front] Solanki|[,0]) THEN writescreen '?Error in file: Solanki-Tapping.ndf' ; writescreen 'Solanki-Tapping data is provided for an irregular date series - there must be constant years difference between data in the series' ; input_error := l ; ELSE Solanki_delta_t := first Solanki_delta_t ; ENDIF ; interpolation_setup ; % I can't think of any checks after initial program debugging ; % Complete the input data by doing extra calculations ; civilisation_period := -11000 + (1000 * tell 12) ; NOT input_error } # Insolation.ndf 17Feb07 www.BillHowell.ca Warning: Use Courier constant-width font and no-line-wrap to view this file! Otherwise the code doesn't line up and it's hard to read. # Original code from: J. Laskar, P. Robutel, F. Joutel, M. Gastineau, A.C.M. Correia, B. Levrard 2004 "A long term numerical solution for the insolation quantities of the Earth" Astronomy & Astrophysics manuscript no. La 2004 May 23, 2004 http://www.imcce.fr/fr/presentation/equipes/ASD/preprints/prep.2004/La_2004_prep.pdf # Laskar, J., Joutel, F., Boudin, F. 1993 "Orbital, precessional and insolation quantities for the Earth from -20 Myr to + 10Myr" Astron. Astrophys. 270, 522 http://www.imcce.fr/Equipes/ASD/insola/earth/La93/README.TXT #**************************************** # Laskar etal - insola.f ported to the Q'Nial programming language by www.BillHowell.ca * Howell: View this file with constant-width font (eg Courier) and word wrap turned off * Much of the original content and references is retained. * * * Computation of the various insolation functions * * this program uses * insolsub.f : subroutines for the computation * of insolation * prepinsol.f : subroutines for the binary temporary file * prepsub.f : subroutines for computing elements at * specified time * insola.par : parameters for insola.f * nomclimaneg , * nomclimapos : ASCII files supplied * which contains the necessary * orbital and precession quantities * * * * the following parameters are taken from the NAMELIST insola.par * nomclimapos : ASCII file for the elements t,e,eps, psibar on * positive time. * nomclimaneg : ASCII file for the elements t,e,eps, psibar on * negative time. * nominsolbin : temporary binary file for the elements t,e,eps, * cos(psibar), sin(psibar) * (variables interne = nomfich) * pas : stepsize (years) * datefin : final time (million of years) * datedebut : starting time (million of years) * so : solar constant W*m^-2 * * other parameters: * * nbf : 4 (number of variables) * nbel1 : 4 + 1 * nacd : numbers of n-uplets in a record * * (c) Astronomie et Systemes Dynamiques/ institut de Mecanique Celeste (2003) * J. Laskar, F. Joutel, M. Gastineau * * version 0.82 (7/7/1993) * version 0.83 (27/1/93) * version la2001b4 (26/06/2001) Mickael : changement du fichier de parametres * version la2001b9 (12/02/2002) Mickael : ajout dans le fichier de parametres * version la2003 (04/03/2003) Mickael : fichier de parametres * version la2003 (04/09/2003) Mickael : creation du fichier binaires * (05/11/2003) Mickael : changement du common date en ddate * version la2004 (26/05/2004) Mickael : correction qromb (insolsub.f) #**************************************** #********************************************* * Laskar etal - insolsub.ndf * * subroutines for the computation of insolation * J. Laskar, F. Joutel, F. Boudin * (c) Astronomie et Systemes Dynamiques / Institut de Mecanique Celeste (2003) * * 05 / 11 / 2001 Mickael : modification de l'appel de polint dans qromb * (mauvais type fourni) * 22 / 09 / 2003 Mickael : utilisation de e et pibarh au lieu de ak,ah, psi * 05 / 11 / 2003 Mickael : changement du common date en ddate d * 26 / 05 / 2004 Mickael : correction qromb suivant numerical recipes 2nd ed. * * traduction en Q'Nial de www.qnial.com - Michael, Jenkins & team at Queen's University, Canada * par www.BillHowell.ca février 2007 * #********************************************* #*************************************** * Description of routines, Howell Feb07 * * e - excentricité * eps - l'obliquité * hl - longitude vraie , moyvrai(hlm, e, pibar) * hlm - longitude moyenne, vraimoy(hl , e, pibar) * phi - latitude du point sur la terre * pibar - longitude du perihelie !&!comptee a partir de l'equinoxe de la date + pi (repere geocentrique)!&! * wd - longitude vraie du soleil (en radian) compte a partir de l'equinoxe de la date * * * * vraimoy (hl,e,pibar,hlm) - calcule la longitude moyenne (hlm) * moyvrai (hlm,e,pibar,hl) - calcule la longitude vraie (hl) * cwj (wd,e,pibar,eps,phi) - INSOLATION JOURNALIERE D'UN POINT DE LATITUDE DONNEE * wmcal (mois,e,eps,pibarh,phi) - INSOLATION MENSUELLE D'UN POINT DE LATITUDE DONNEE * wjour (date,e,eps,pibarh,phi) - INSOLATION JOURNALIERE D'UN POINT DE LATITUDE DONNEE * wjcal (datecal,e,eps,pibarh,phi) - INSOLATION JOURNALIERE D'UN POINT DE LATITUDE DONNEE * Telinsol (tps) - tableau des éléments orbitaux (k,h,eps,phi) pour la date tps * (Howell: telinsol not called in Laskar program !!&!!!&!!) NOT called by prepinsol.f nor prepsub.f !!!!! * Global_annual_avg_insolation (e) - INSOLATION ANNUELLE MOYENNEE * * * From standard libraries at !&!l'Observatoire de Paris!&! * (c) Astronomie et Systemes Dynamiques/ institut de Mecanique Celeste (2003) * * QROMB (FUNC,A,B,SS) - ROMBERG INTEGRATION * POLINT (XA, YA, N, X) - used by qromb * TRAPZD (FUNC A B N) - used by qromb * * ******************************************* #********************************* # globals defined here, or used in this file only ; e_du := 0.0 ; e_du_v := f3 ; eps_du := 0.0 ; eps_du_v := f3 ; nbfmax := 20 ; nbfmax_v := f3 ; q_dy := 0.0 ; q_dy_v := f3 ; % for rhomberg routines ; q_nmax := 10 ; q_nmax_v := f3 ; % for rhomberg routines ; phi_du := 0.0 ; phi_du_v := f3 ; pibar_du := 0.0 ; pibar_du_v := f3 ; % for rhomberg integration routines ; q_c := q_nmax reshape 0.0 ; q_c_v := f3 ; % for rhomberg routines ; q_d := q_nmax reshape 0.0 ; q_d_v := f3 ; % for rhomberg routines ; q_xa := q_nmax reshape 0.0 ; q_xa_v := f3 ; % for rhomberg routines ; q_ya := nbfmax q_nmax reshape 0.0 ; q_ya_v := f3 ; % for rhomberg routines ; cwj_o := f3 ; civilization_insolation_o := f3 ; Global_annual_avg_insolation_o := f3 ; Isolation_mensuel_o := f3 ; moyvrai_o := f3 ; QROMB_o := f3 ; POLINT_o := f3 ; Telinsol_o := f3 ; TRAPZD_o := f3 ; vraimoy_o := f3 ; wmcal_o := f3 ; wjour_o := f3 ; wjcal_o := f3 ; # EXTERNAL definitions to allow this file to load alone civilisation_period IS external variable ; civilisation_period_x := f3 ; fichout IS external variable ; fichout_x := f3 ; intrp_tseries IS external variable ; intrp_tseries_x := f3 ; Milankovic IS external variable ; Milankovic_x := f3 ; mois IS external variable ; mois_x := f3 ; n_intrpsteps IS external variable ; n_intrpsteps_x := f3 ; phi IS external variable ; phi_x := f3 ; so IS external variable ; so_x := f3 ; interpolation_for_interval IS external operation ; interpolation_for_interval_x := f3 ; # ********************************************* c subroutine vraimoy(hl,e,pibar,hlm) c c calcule la longitude moyenne (hlm) en fonction de la longitude c vraie (hl), de l'excentricite (e) et de la longitude du perihelie c (pibar) c (c) ASD / BDL (1 / 10 / 92) c correction 27 / 1 / 93 # ********************************************* vraimoy IS OP hl e pibar { beta := (1 - (e * e)) power 0.5 ; eM := hl - pibar ; hl - (2 * (e / 2 + (e power 3 / 8) ) * (1 + beta) * sin eM ) - ( (e * e / 4 ) * (0.5 + beta) * sin (2 * eM) ) + ( (e power 3 ) * (1 / 3 + beta) * sin (3 * eM) ) } # ********************************************* c subroutine moyvrai(hlm,e,pibar,hl) c c calcule la longitude vraie (hl) en fonction de la longitude c moyenne (hlm), de l'excentricite (e) et de la longitude du c perihelie (pibar) c (c) ASD / BDL (1 / 10 / 92) # ********************************************* moyvrai IS OP hlm e pibar { eM := hlm - pibar ; hlm + ( ( 2 * e - (e power 3 / 4)) * sin eM ) + ( 5 / 4 * (e * e) * sin (2 * eM ) ) + ( 13 / 12 * (e power 3) * sin (3 * eM ) ) } # ********************************************* c SUBROUTINE cwj(wd,e,pibar,eps,phi,w) c c INSOLATION JOURNALIERE D'UN POINT DE LATITUDE DONNEE c c A L'ENTREE: c wd : la longitude vraie du soleil (en radian) compte a partir c de l'equinoxe de la date c e : excentricite c pibar: longitude du perihelie comptee a partir de l'equinoxe de la c date + pi (repere geocentrique) c eps: l'obliquite c phi: latitude du point sur la terre c c EN SORTIE: c w:l'insolation journaliere d'un point de latitude donnee c (c) ASD / BDL (1 / 10 / 92) # ********************************************* cwj IS OP wd e pibar eps phi { v := wd - pibar ; rho := (1 - (e * e)) / (1 + (e * cos v) ) ; % Distance terre - soleil (rho): ; delta := arcsin (sin eps * sin wd) ; % Declinaison du soleil (delta): ; aux := pi / 2 - abs delta ; IF (( - aux < phi) and (phi < aux)) THEN % LATITUDES DE LEVER ET COUCHER DE SOLEIL: ; cho := - (tan phi) * (tan delta) ; ho := arccos cho ; % Angle horaire du lever et coucher du soleil (ho): ; ( ( ho * (sin phi) * (sin delta) ) + ( (cos phi) * (cos delta) * (sin ho) ) ) * so / (pi * (rho * rho) ) ELSE a1 := pi / 2 - delta ; a2 := pi / 2 + delta ; IF ((phi >= a1) or (phi <= ( - a2)) ) THEN % LATITUDES SANS COUCHER: ; so * (sin phi) * (sin delta) / (rho * rho) ELSE % LATITUDES SANS LEVER: ; IF ((phi <= ( - a1)) OR (phi >= a2) ) THEN 0 ENDIF ENDIF ENDIF } # ********************************************* * Isolation_mensuel # ********************************************* * Isolation_mensuel IS OP hlm { cwj (moyvrai hlm e_du pibar_du) e_du pibar_du eps_du phi_du } # ********************************************* * c SUBROUTINE wjour(date,e,eps,pibarh,phi,w) c c INSOLATION JOURNALIERE D'UN POINT DE LATITUDE DONNEE c c A L'ENTREE: c date : la longitude vraie du soleil (en degres) comptee a partir c de l'equinoxe vrai c e,eps,pibarh: les elements orbitaux pour la date tps c e : excentricite c eps:l'obliquite c pibarh: longitude du perihelie comptee a partir de l'equinoxe de la c date c phi: latitude du point sur la terre c c c ---------------------------------------------------------------------- - c CALCUL D'INSOLATION c ---------------------------------------------------------------------- - c c longitude du perihelie par rapport a l'equinoxe de reference(pitild), c par rapport a l'equinoxe de la date(pibar), excentricite(e): c c c Longitude vraie comptee a partir de l'equinoxe vrai et c anomalie vraie (wd,v): c c c EN SORTIE: c w:l'insolation journaliere d'un point de latitude donnee c (c) ASD / IMC (22 / 09 / 2003) # ********************************************* * wjour IS OP date e eps pibarh phi { pibar := pibarh + pi ; hl := date * pi / 180 ; cwj hl e pibar eps phi } # ********************************************* * c SUBROUTINE wjcal(datecal,e,eps,pibarh,phi,w) c c INSOLATION JOURNALIERE D'UN POINT DE LATITUDE DONNEE c c A L'ENTREE: c datecal: la longitude moyenne (en degres) du soleil comptee a partir c de l'equinoxe de la date c e,eps,pibarh: les elements orbitaux pour la date tps c e : excentricite c eps:l'obliquite c pibarh: longitude du perihelie comptee a partir de l'equinoxe de la c date c phi: latitude du point sur la terre c c c ---------------------------------------------------------------------- - c CALCUL D'INSOLATION c ---------------------------------------------------------------------- - c c longitude du perihelie par rapport a l'equinoxe de reference(pitild), c par rapport a l'equinoxe de la date(pibar), excentricite(e): c c c EN SORTIE: c w:l'insolation journaliere d'un point de latitude donnee c (c) ASD / BDL (1 / 10 / 92) # ********************************************* * wjcal IS OP datecal e eps pibarh phi { pibar := pibarh + pi ; hlm := datecal * pi / 180 ; % calcul de la longitude vraie ; hlm0 := vraimoy 0. e pibar ; % longitude moyenne au 21 mars ; hlm := hlm0 + hlm ; % longitude moyenne a la date ; wd := moyvrai hlm e pibar ; % longitude vraie a la date ; cwj wd e pibar eps phi } # ********************************************* c SUBROUTINE Global_annual_avg_insolation(e,w) c c INSOLATION ANNUELLE MOYENNEE c c A L"ENTREE: c e:excentricite pour une date donnee c c ---------------------------------------------------------------------- - c CALCUL D'INSOLATION (w) c ---------------------------------------------------------------------- - c Carre de la distance terre - soleil moyennee sur une revolution (rho): c c c EN SORTIE: c w:insolation annuelle moyenne c (c) ASD / BDL (1 / 10 / 92) # ********************************************* Global_annual_avg_insolation IS OP e { so / (4 * sqrt(1 - e power 2)) } # ************************************************* c SUBROUTINE QROMB(FUNC,A,B,SS) c c NUMERICAL RECIPES c ROMBERG INTEGRATION # ************************************************* c 05 / 11 / 2001 Mickael : ajout de 0 (au lieu de 0.) pour l'appel de polint c 26 / 05 / 2004 Mickael : correction du test IF (ABS(DSS).LT. en .LE. c correction effectuee dans "numerical recipes, second edition" POLINT IS OP N X { NONLOCAL q_c q_d q_dy q_xa q_y q_ya q_nmax ; % ; NS := 1 ; DIF := ABS(X - q_xa@1) ; FOR I WITH count N DO DIFT := ABS (X - q_xa@I) ; IF DIFT < DIF THEN NS := I ; DIF := DIFT ; ENDIF ; q_C@I := q_ya@I ; q_D@I := q_ya@I ; ENDFOR ; q_y := q_ya@NS ; NS := NS - 1 ; FOR M WITH count (N - 1) DO FOR I WITH count (N - M) DO HO := q_xa@I - X ; HP := q_xa@[I + M] - X ; W := q_c@[I + 1] - q_d@I ; DEN := HO - HP ; IF DEN = 0.0 THEN writescreen ' pause polint - den=0 error where interpolated point is at known point' ; ENDIF ; % what does this do!&! ; DEN := W / DEN ; q_d@I := HP * DEN ; q_c@I := HO * DEN ; ENDFOR ; IF ( (2 * NS) < (N - M) ) THEN q_dy := q_c@[NS + 1] ; ELSE q_dy := q_d@NS ; NS := NS - 1 ; ENDIF ; q_y:= q_y+ q_dy ; ENDFOR ; % return 2 results ; q_y q_dy } TRAPZD IS OP FUNC A B N { % SAVE IT ; % what does IT do!&! intermediate save in file!&! ; IF N = 1 THEN S := 0.5 * (B - A) * ((apply FUNC A) + (apply FUNC B)) ; IT := 1 ; ELSE TNM := IT ; DEL := (B - A) / TNM ; X := A + (0.5 * DEL) ; SUMT := 0. ; FOR J WITH count IT DO SUMT := SUMT + (FUNC X) ; X := X + DEL ; ENDFOR ; S := 0.5 * (S + ( (B - A) * SUMT / TNM) ); IT := 2 * IT ; ENDIF ; S } QROMB IS OP FUNC A B { error_qromb := o ; EPS := 1.e-6 ; JMAX := 20 ; JMAXP := JMAX + 1 ; K := 5 ; KM := 4 ; S := JMAXP reshape 0.0 ; H := S ; H@1 := 1. ; FOR J WITH count JMAX DO S@J := TRAPZD FUNC A B J ; IF J >= K THEN Lx := J - KM ; SS DSS := POLINT H@Lx S@Lx K 0.0 ; IF (ABS DSS <= (EPS * ABS SS) ) THEN error_qromb := l ; ENDIF ; ENDIF ; S@(J + 1) := S@J ; H@(J + 1) := 0.25 * H@J ; ENDFOR ; IF error_qromb = l THEN Name := Readscreen 'QROMB: Too many steps.' ELSE SS ENDIF } # ********************************************* * c SUBROUTINE wmcal(mois,e,eps,pibarh,phi,w) c c INSOLATION MENSUELLE D'UN POINT DE LATITUDE DONNEE c c A L'ENTREE: c mois : le numero du mois (l'annee est divisee en 12 mois de 30 c degres et le numero 3 correspond a fin fevrier - fin mars) c e,eps,pibarh: les elements orbitaux pour la date tps c e : excentricite c eps:l'obliquite c pibarh: longitude du perihelie comptee a partir de l'equinoxe de la c date c phi: latitude du point sur la terre c c c ---------------------------------------------------------------------- - c CALCUL D'INSOLATION c ---------------------------------------------------------------------- - c c longitude du perihelie par rapport a l'equinoxe de reference(pitild), c par rapport a l'equinoxe de la date(pitbar), excentricite(e): c c Longitude moyenne au 21 mars c c c EN SORTIE: c w:l'insolation mensuelle d'un point de latitude donnee c (c) ASD / IMC (22 / 09 / 2003) # ********************************************* * # Howell - this probably isn't useful (including QROMB) forthe very short timeframes that I am looking at wmcal IS OP mois e eps pibarh phi { NONLOCAL phi_du eps_du e_du pibar_du ; % keep these intermediates, as I don't know if conflicts will arise by making e eps puibarh etc global ; phi_du := phi ; eps_du := eps ; e_du := e ; % globals so that QROMB & Isolation_mensuel can work ; pibar_du := pibarh + pi ; % ; hlm0 := vraimoy 0 e pibar_du ; % longitude moyenne (hlm0) au debut de l'annee!&! ; hlm1 := hlm0 + (mois - 4) * pi * 30 / 180 ; % longitude moyenne (hlm1) au debut du mois ; hlm2 := hlm1 + pi * 30 / 180 ; % longitude moyenne (hlm2) au fin du mois ; % calcul de l'insolation moyenne a l'aide de la methode de Romdberg ; (qromb "Isolation_mensuel hlm1 hlm2) / 30 / pi * 180 } civilization_insolation IS OP fichout latitude mois { LOCAL interval_inputs interval_outpts ; % ; % Solanki ; yr := 0 ; spot := 1 ; irr := 2 ; STf := 3 ; % Milankovic ; yr := 0 ; ecc := 1 ; obl := 2 ; cosp := 3 ; sinp := 4 ; prec := 5 ; % interval_parms ; yr := 0 ; ecc := 1 ; obl := 2 ; cosp := 3 ; sinp := 4 ; prec := 5 ; % interval_inputs ; mois := 0 ; ecc := 1 ; obl := 2 ; pr := 3 ; lat := 4 ; % interval_outpts ; yr := 0 ; ecc := 1 ; obl := 2 ; cosp := 3 ; sinp := 4 ; prec := 5 ; irr := 6 ; inraw := 7 ; insol := 8 ; % ; interval_inputs := n_intrpsteps 5 reshape 0.0 ; % initial assignments same for full problem set ; interval_inputs|[,mois] := n_intrpsteps reshape mois ; % mean monthly, eg 7 => 21 june - 20 july ; interval_inputs|[,lat ] := n_intrpsteps reshape lat ; % latitude du point sur la terre ; interval_outpts := n_intrpsteps 9 reshape 0.0 ; % ; fichout_number := open fichout "w ; fichout_number EACHRIGHT writefile ' ' fichout (link 'Averages for the month of: ' string mois) '' 'Year eccentricity obliquity precession cos_p sin_p irradiance insol_raw insolation' ; % ; FOR interval_start_year WITH civilisation_period DO interval_years := interval_start_year + intrp_tseries ; interval_outpts|[, yr ecc obl cosp sinp] := interpolation_for_interval interval_start_year ; % ; ST_start_index := find interval_start_year Solanki|[,yr] ; interval_STfs := Solanki|[ST_start_index + intrp_tseries, STf] ; % ; interval_inputs|[,ecc obl] := interval_outpts|[,ecc obl] ; interval_inputs|[,pr] := EACH arctan (OUTER / interval_outpts|[,sinp] interval_outpts|[,cosp]) ; interval_outpts|[,inraw] := EACH wmcal rows interval_inputs ; interval_outpts|[,insol] := OUTER * interval_outpts|[,inraw] interval_STfs ; % ; EACH writefile EACH string rows interval_outpts ; ENDFOR ; close fichout_number ; } ?undefined identifier: FIND INTERVAL_START_YEAR SOLANKI <***> | [ , # enddoc errors found: 1 # Interpolation.ndf 17Feb07 www.BillHowell.ca port from Fortran to QNial Warning: Use Courier constant-width font and no-line-wrap to view this file! Otherwise the code doesn't line up and it's hard to read. # Original code from: J. Laskar, P. Robutel, F. Joutel, M. Gastineau, A.C.M. Correia, B. Levrard 2004 "A long term numerical solution for the insolation quantities of the Earth" Astronomy & Astrophysics manuscript no. La 2004 May 23, 2004 http://www.imcce.fr/fr/presentation/equipes/ASD/preprints/prep.2004/La_2004_prep.pdf # Laskar, J., Joutel, F., Boudin, F. 1993 "Orbital, precessional and insolation quantities for the Earth from -20 Myr to + 10Myr" Astron. Astrophys. 270, 522 http://www.imcce.fr/Equipes/ASD/insola/earth/La93/README.TXT #********************************* Introduction # This set of routines has been simplified from Laskar etal and adapted more for the Holocene period of geology (roughly the post glaciation period leading into civilization and up to the present day). # Because the ultimate intent is to incorporate time series for solar, geomagnetic and galactic drivers of climate and disease (and many other posiible effects, like war, agriculture, the economy, arts, sciences etc) a time resolution on the order of one year is desired. However, Solanki's proxy time seires for sunspots, which was constructed from geologic data including 10Be and 14C, only has a 5 year time resolution at present, so I'll use that time resolution as a fixed parameter (in a way to make it easy for someone else to change it). # Solanki, S.K., et al. 2005. 11,000 Year Sunspot Number Reconstruction. IGBP PAGES/World Data Center for Paleoclimatology Data Contribution Series #2005-015. NOAA/NGDC Paleoclimatology Program, Boulder CO, USA. # As the Laskar etal Milankovic parameters (eccentricity, obliquity, precession, but not Earth's bobbing in and out of the plane of its orbit, etc etc) are provided every 1,000 years, there is a constant 1,000 / 5 = 200 data points of interpolation between each provided Milankovic point. We also want a match between "current dates" for the Laskar and Solanki data series - but that will have to come later. # It is my impression that the Laskar interpolation program allowed for Milankovic series with irregular time intervals. But the provided Milankovic time series is for a relular 1,000 year time interval, so for this assumption we can simplify the calcualtions somewhat. # Obviously, there is a bit of a "drift" in using the Solanki and Laskar time series "before present", as opposed to using a fixed reference date (eg 2000 AD, which somehow seems politically incorrect to say, such is the nonsense we live with). # # #************************************************* temporary - for helping to write routines... f2 := 'Interpolation.ndf' ; #********************************* # globals defined here, or used in this file only ; Holocene_start := -25000 ; Holocene_start_v := f2 ; % This is simply the start ofthe Milankovic series for interpolation, not the geological period!! ; Holocene_end := 2000 ; Holocene_end_v := f2 ; % 2,000 years into the future ; intrp_factors_v := f2 ; intrp_tdata := 0.0 ; intrp_tdata_v := f2 ; intrp_tseries := 0.0 ; intrp_tseries_v := f2 ; n_intrpsteps := 0.0 ; n_intrpsteps_v := f2 ; n_poly := 8 ; n_poly_v := f2 ; % 8th order polynomial fit ; nbfmax := 20 ; nbfmax_v := f2 ; nmax := 100 ; nmax_v := f2 ; n_vars_e := f2 ; % # of variables -= 4 (eccentricity, obliquity, cos precess, sin precess) ; facts_pre&aft_interval := 0.0 ; facts_pre&aft_interval_v := f2 ; ps_c := nmax reshape 0.0 ; ps_c_v := f2 ; ps_d := nmax reshape 0.0 ; ps_d_v := f2 ; ps_dy := nbfmax reshape 0.0 ; ps_dy_v := f2 ; interpolation_for_interval_o := f2 ; interpolation_setup_o := f2 ; pint_o := f2 ; # EXTERNAL definitions to allow this file to load alone, some defined in header file Milankovic IS external variable ; Milankovic_x := f2 ; Milankovic_rows IS external variable ; Milankovic_rows_x := f2 ; Milankovic_cols IS external variable ; Milankovic_cols_x := f2 ; Milankovic_delta_t IS external variable ; Milankovic_delta_t_x := f2 ; Solanki_delta_t IS external variable ; Solanki_delta_t_x := f2 ; #*************************************************** % SUBROUTINE PINT(xa,ya,n,x,nbf,y) ; % % % A partir de fonctions tabulees (xa,ya),ce sous - programme construit ; % un polynome de degre (n - 1) permettant de calculer la valeur de la ; % fonction, y, au point, x. ; % (cf algorithme de Neville, Numerical recipes) ; % ; % A L'ENTREE: ; % xa(n):tableau de dates ; % ya(nbf,n):valeurs des fonctions aux dates xa(n) ; % n:nombre de points necessaires a l'interpolation ; % x:date pour laquelle on effectue l'interpolation ; % nbf:nombre de fonctions a interpoler ; % ; % EN SORTIE: ; % y(nbf):tableau des valeurs interpolees a la date x ; % ; % (c) ASD / IMC (2001) ; pint IS OP x { % Milankovic ; yr := 0 ; ecc := 1 ; obl := 2 ; cosp := 3 ; sinp := 4 ; prec := 5 ; % ; ps_xa := facts_pre&aft_interval|[,yr ] ; ps_ya := facts_pre&aft_interval|[,ecc obl cosp sinp ] ; ps_y := shape ps_ya reshape 0.0 ; n := n_poly ; nbf := n_poly ; % ; FOR nvar WITH tell nbf DO ns := 1 ; dif := abs (x - ps_xa@1) ; FOR i WITH tell n DO dift := abs (x - ps_xa@i) ; IF (dift < dif) THEN ns := i ; dif := dift ; ENDIF ; ps_c@i := ps_ya@[nvar,i] ; ps_d@i := ps_ya@[nvar,i] ; ENDFOR ; ps_y@nvar := ps_ya@[nvar,ns] ; ns := ns - 1 ; FOR m WITH tell (n - 1) DO FOR i WITH tell (n - m) DO ho := ps_xa@i - x ; hp := ps_xa@[i + m] - x ; w := ps_c@[i + 1] - ps_d@i ; den := ho - hp ; IF (den = 0) THEN writescreen ' pause pint - den=0 error where interpolated point is at known point' ; break ; ENDIF ; den := w / den ; ps_d@i := hp * den ; ps_c@i := ho * den ; ENDFOR ; IF (2 * ns) < (n - m) THEN ps_dy@nvar := ps_c@[ns + 1] ; ELSE ps_dy@nvar := ps_d@ns ; ns := ns - 1 ; ENDIF ; ps_y@nvar := ps_y@nvar + ps_dy@nvar ; ENDFOR ; ENDFOR ; ps_y } #*************************************************** interpolation_setup IS { n_intrpsteps := floor ( Milankovic_delta_t / Solanki_delta_t ) ; % includes start of interval - year of given Milankovic data (must check result consistency!) ; intrp_tseries := Solanki_delta_t * tell n_intrpsteps ; intrp_tdata := OUTER + ( intrp_tseries ) ( (- (tell n_poly) 4 ) * Milankovic_delta_t ) ; } #*************************************************** Returns a table of interpolated Milankovic_factors_for_interval|[,0] - year eccentricity obliquity cos(precession) sin(precession) At some point verify that the years match (data vs interpolation for 1st year)!!! interpolation_for_interval IS OP interval_start_year { % Milankovic columns ; yr := 0 ; ecc := 1 ; obl := 2 ; cosp := 3 ; sinp := 4 ; prec := 5 ; % ; startYear_Index := find interval_start_year Milankovic|[,yr] ; startYear_data := Milankovic|[floor (startYear_Index + (n_poly / 2) - 1),] ; facts_pre&aft_interval := Milankovic|[startYear_Index + tell n_poly, yr ecc obl cosp sinp] ; EACH pint intrp_tseries } # enddoc # Milankovic data - estimates of the Milankovic parameters for the period of civilizations Civilization Milankovic precess, obliq, eccent.txt 11Feb07 www.BillHowell.ca Warning: Use Courier constant-width font and no-line-wrap to view this file! Otherwise the code doesn't line up and it's hard to read. # Original code from: J. Laskar, P. Robutel, F. Joutel, M. Gastineau, A.C.M. Correia, B. Levrard 2004 "A long term numerical solution for the insolation quantities of the Earth" Astronomy & Astrophysics manuscript no. La 2004 May 23, 2004 http://www.imcce.fr/fr/presentation/equipes/ASD/preprints/prep.2004/La_2004_prep.pdf # Laskar, J., Joutel, F., Boudin, F. 1993 "Orbital, precessional and insolation quantities for the Earth from -20 Myr to + 10Myr" Astron. Astrophys. 270, 522 http://www.imcce.fr/Equipes/ASD/insola/earth/La93/README.TXT #********************************* # globals defined here, or used in this file only ; date_conversion := 1000 ; date_conversion_v := f10 ; % convert dates to years (eg from My or ky) ; Milankovic := 0.0 ; Milankovic_v := f10 ; % input below ; Milankovic_cols := 6 ; Milankovic_cols_v := f10 ; % important for reading data!! ; Milankovic_rows := 0.0 ; Milankovic_rows_v := f10 ; % calc'd below ; Milankovic_titles := ' ' ; Milankovic_titles_v := f10 ; n_vars_v := f10 ; so := 1368 ; so_v := f10 ; % solar constant (W*m^-2), this can be modified in the input file! ; % there aren't any operators defined in this file ; # EXTERNAL definitions to allow this file to load alone # none required... # Some nomenclature used by Laskar etal e or k % eccentricity ; eps or h % obliquity ; pibar % precession - !&!l'angle entre longitude de l'equinox et perihelion!&! ; q % cos(precession) ; p % cos(precession) ; pas % stepsize (years) for interpolations, must be greater than 1 year ; % defined in 'Civilization problem setup.ndf' ; #******************************************************** # Solar insolation parameters from Laskar etal # date (kyears) - relative to current time. Note that this is translated into years for later calculations but it's easiest to retain the ky inpout style for the user for the Holocene period - Milankovic dates must start with earliest and end with last!! (time goes forward) - NOTE:, I'm not sure exactly when "present time" is!!! assume ~2005 (like SOlaki basis I'm using) # eccentricity (!&!units!&!) - measure of degree to which Earth's orbit around the sun is non-circular - 90, 100, 110 ky periodicity (taken to be ~100 ky) # obliquity (radians) - angle of Earth's axis with respect to its plane of orbit - ~40 ky periodicity # precession (radians) - direction of Earth's axis tilt with respect to !&!today's orientation versus the center of the universe!&! - ~20 ky periodicity # cos(precession) - all zeros, waiting for the results to be put in # sin(precession) - all zeros, waiting for the results to be put in # "Invisible column" - a tab and a space finish each row - otherwise the insolation & following year are read as one!! # initially the precession is in terms of an angle (radians), then it is converted to sin & cos terms. Also from kYears to Years. Oops, when I copied overthe data it went from 16 to 9 digit accuracy (double to single precision) I'll have to fix later. # solar constant (W*m^-2), already defined in external.ndf, but this can be modified below! The next aim is to use Solanki's solar irradiance values for the Holocene. # copy this into operators using Milankovic % Milankovic columns ; yr_col := 0 ; ecc_col := 1 ; obl_col := 2 ; cosp_col := 3 ; sinp_col := 4 ; prec_col := 5 ; Milankovic_titles := 'YearsBP' 'eccentricity' 'obliquity' 'cos(precession)' 'sin(precession)' 'precession' ; Milankovic := -25 0.0178680000 0.3911080000 0.0 0.0 0.8712670000 -24 0.0180950000 0.3931010000 0.0 0.0 1.1553320000 -23 0.0183410000 0.3954330000 0.0 0.0 1.4405280000 -22 0.0186170000 0.3980270000 0.0 0.0 1.7259070000 -21 0.0188350000 0.4008000000 0.0 0.0 2.0112180000 -20 0.0190030000 0.4036770000 0.0 0.0 2.2983900000 -19 0.0192230000 0.4065770000 0.0 0.0 2.5814010000 -18 0.0193470000 0.4094130000 0.0 0.0 2.8676430000 -17 0.0194870000 0.4121230000 0.0 0.0 3.1539570000 -16 0.0195600000 0.4146340000 0.0 0.0 3.4367300000 -15 0.0196350000 0.4168870000 0.0 0.0 3.7227540000 -14 0.0196790000 0.4188350000 0.0 0.0 4.0084460000 -13 0.0196510000 0.4204360000 0.0 0.0 4.2923470000 -12 0.0196130000 0.4216580000 0.0 0.0 4.5794320000 -11 0.0195970000 0.4224860000 0.0 0.0 4.8634580000 -10 0.0194250000 0.4229060000 0.0 0.0 5.1533730000 -9 0.0193250000 0.4229180000 0.0 0.0 5.4383550000 -8 0.0191800000 0.4225360000 0.0 0.0 5.7272450000 -7 0.0189370000 0.4217770000 0.0 0.0 6.0186410000 -6 0.0187480000 0.4206670000 0.0 0.0 0.0245740000 -5 0.0184530000 0.4192470000 0.0 0.0 0.3155370000 -4 0.0182110000 0.4175550000 0.0 0.0 0.6113810000 -3 0.0178460000 0.4156390000 0.0 0.0 0.9034940000 -2 0.0174970000 0.4135550000 0.0 0.0 1.2000240000 -1 0.0171610000 0.4113530000 0.0 0.0 1.4979280000 0 0.0167020000 0.4090930000 0.0 0.0 1.7962570000 1 0.0162751900 0.4068337800 0.0 0.0 2.0985607000 2 0.0158501400 0.4046314300 0.0 0.0 2.3976566000 3 0.0153397000 0.4025410200 0.0 0.0 2.7053148000 4 0.0148354800 0.4006169700 0.0 0.0 3.0096091000 5 0.0143164700 0.3989009600 0.0 0.0 3.3134127000 6 0.0137745200 0.3974365500 0.0 0.0 3.6261098000 7 0.0132086400 0.3962542200 0.0 0.0 3.9336626000 8 0.0125777100 0.3953768400 0.0 0.0 4.2441453000 9 0.0120396300 0.3948181800 0.0 0.0 4.5577560000 10 0.0113841500 0.3945853300 0.0 0.0 4.8728309000 11 0.0107391500 0.3946720400 0.0 0.0 5.1874266000 12 0.0101621900 0.3950691200 0.0 0.0 5.5067319000 13 0.0094765790 0.3957572900 0.0 0.0 5.8251114000 14 0.0088686560 0.3967117900 0.0 0.0 6.1514530000 15 0.0082116980 0.3979020200 0.0 0.0 0.1866843700 ; Milankovic_rows := first floor ((shape Milankovic) / Milankovic_cols ) ; Milankovic := Milankovic_rows Milankovic_cols reshape Milankovic ; n_vars := (second shape Milankovic) - 1 ; % later merge years into this ; # should double-check to make sure "Milankovic_rows Milankovic_cols" have been changed if data changes # enddoc # Solanki-Tapping.ndf Sunspot estimates by Solanki (10Be and 14C based) Irradiance model of Tapping #********************************************************************* 11,000 Year Sunspot Number Reconstruction. "Solanki, S.K., et al. 2005. 11,000 Year Sunspot Number Reconstruction. IGBP PAGES/World Data Center for Paleoclimatology Data Contribution Series #2005-015. NOAA/NGDC Paleoclimatology Program, Boulder CO, USA. " # DESCRIPTION: The series of reconstructed 10-yr averaged sunspot numbers with their 68% uncertainty. Years are given BP (before present), i.e. the calendar AD year, Yad, is related to the BP year, Ybp, as Yc=1950-Ybp. The tabulated years correspond to centers of the corresponding 10-year intervals. Negative values are artifacts and are consistent with zero within the error limits. # DATA: Column 1: Year Before Present (1950 AD) Column 2: Reconstructed Sunspot Number Column 3: Sigma error (68% uncertainty) # #********************************* Howell - # Here YearAD has a +105 shift to make "present time" = 2005 AD (i.e. forst data in 1895) (just to make it easy for me to align Solanki and Laskar datasets (sheesh, I'm lazy) # The intent is to provide Tapping's model of Total irradiance as a function of sunspots, which will come from Solanki's time series. # The Solanki_Tapping time series MUST cover the same "interval" as the "civilisation period" used for the Lasskar insolation calculations although extra years before or after are OK. Later improvements may relax that constraint. Obviously, data since 1895 for sunspots would be helpful!!! # Both the Solanki_Tapping series and "Mialnkovic_delta_t" must have the same interval time step. Initially it has been taken to be 10 years, as with the Solanki data series. # Years that don't have Solanki sunspot data have simply been provided with the "solar constant" = 1368 Wm^-2 #********************************* # globals defined here, or used in this file only ; f11 := 'Solanki-Tapping.ndf' ; Solanki := 0.0 ; Solanki_v := f11 ; % input below ; Solanki_cols := 4 ; Solanki_cols_v := f11 ; % important for reading data!! ; Solanki_rows := 0.0 ; Solanki_rows_v := f11 ; % calc'd below ; Solanki_titles := 'YearBP_shifted' 'Sunspot #' ; Solanki_titles_v := f11 ; Solanki_Tapping_irradiance := 0.0 ; % input below (or will be when I get this all working) ; Solanki_Tapping_irradiance_v := f11 ; % there aren't any operators defined in this file ; # EXTERNAL definitions to allow this file to load alone # none required... #********************************* # copy this into start of operatiuons, expressions % for Solanki-Tapping ; y_col := 0 ; spot_col := 1 ; irr_col := 2 ; STf_col := 3 : Solanki_titles := 'YearBP_shifted' 'Solanki sunspot#' 'Tapping_irradiance' 'S-T factor' +--------------+----------------+------------------+----------+ |YearBP_shifted|Solanki sunspot#|Tapping_irradiance|S-T factor| +--------------+----------------+------------------+----------+ Solanki := -11460 37.40 1368.00 1.0000000 -11450 57.10 1368.00 1.0000000 -11440 89.00 1368.00 1.0000000 -11430 92.50 1368.00 1.0000000 -11420 75.20 1368.00 1.0000000 -11410 60.60 1368.00 1.0000000 -11400 43.30 1368.00 1.0000000 -11390 37.70 1368.00 1.0000000 -11380 57.10 1368.00 1.0000000 -11370 71.80 1368.00 1.0000000 -11360 51.90 1368.00 1.0000000 -11350 42.30 1368.00 1.0000000 -11340 54.80 1368.00 1.0000000 -11330 65.00 1368.00 1.0000000 -11320 50.30 1368.00 1.0000000 -11310 38.10 1368.00 1.0000000 -11300 36.40 1368.00 1.0000000 -11290 28.40 1368.00 1.0000000 -11280 23.40 1368.00 1.0000000 -11270 25.00 1368.00 1.0000000 -11260 20.30 1368.00 1.0000000 -11250 14.50 1368.00 1.0000000 -11240 12.30 1368.00 1.0000000 -11230 7.20 1368.00 1.0000000 -11220 3.20 1368.00 1.0000000 -11210 6.50 1368.00 1.0000000 -11200 12.60 1368.00 1.0000000 -11190 15.80 1368.00 1.0000000 -11180 16.20 1368.00 1.0000000 -11170 16.10 1368.00 1.0000000 -11160 16.90 1368.00 1.0000000 -11150 16.60 1368.00 1.0000000 -11140 15.30 1368.00 1.0000000 -11130 12.40 1368.00 1.0000000 -11120 8.50 1368.00 1.0000000 -11110 7.90 1368.00 1.0000000 -11100 8.90 1368.00 1.0000000 -11090 14.00 1368.00 1.0000000 -11080 26.70 1368.00 1.0000000 -11070 32.50 1368.00 1.0000000 -11060 25.80 1368.00 1.0000000 -11050 21.90 1368.00 1.0000000 -11040 23.90 1368.00 1.0000000 -11030 31.80 1368.00 1.0000000 -11020 43.60 1368.00 1.0000000 -11010 57.00 1368.00 1.0000000 -11000 81.10 1368.00 1.0000000 -10990 105.80 1368.00 1.0000000 -10980 93.80 1368.00 1.0000000 -10970 77.60 1368.00 1.0000000 -10960 70.40 1368.00 1.0000000 -10950 67.20 1368.00 1.0000000 -10940 64.00 1368.00 1.0000000 -10930 57.10 1368.00 1.0000000 -10920 49.20 1368.00 1.0000000 -10910 44.50 1368.00 1.0000000 -10900 45.30 1368.00 1.0000000 -10890 58.70 1368.00 1.0000000 -10880 76.20 1368.00 1.0000000 -10870 73.10 1368.00 1.0000000 -10860 64.30 1368.00 1.0000000 -10850 59.50 1368.00 1.0000000 -10840 54.90 1368.00 1.0000000 -10830 53.30 1368.00 1.0000000 -10820 48.60 1368.00 1.0000000 -10810 40.00 1368.00 1.0000000 -10800 32.80 1368.00 1.0000000 -10790 30.70 1368.00 1.0000000 -10780 38.40 1368.00 1.0000000 -10770 49.30 1368.00 1.0000000 -10760 44.30 1368.00 1.0000000 -10750 29.70 1368.00 1.0000000 -10740 23.60 1368.00 1.0000000 -10730 24.40 1368.00 1.0000000 -10720 34.10 1368.00 1.0000000 -10710 43.90 1368.00 1.0000000 -10700 32.90 1368.00 1.0000000 -10690 23.90 1368.00 1.0000000 -10680 36.70 1368.00 1.0000000 -10670 57.60 1368.00 1.0000000 -10660 59.00 1368.00 1.0000000 -10650 45.20 1368.00 1.0000000 -10640 35.10 1368.00 1.0000000 -10630 30.60 1368.00 1.0000000 -10620 27.40 1368.00 1.0000000 -10610 24.20 1368.00 1.0000000 -10600 27.10 1368.00 1.0000000 -10590 34.90 1368.00 1.0000000 -10580 42.20 1368.00 1.0000000 -10570 41.40 1368.00 1.0000000 -10560 32.00 1368.00 1.0000000 -10550 25.40 1368.00 1.0000000 -10540 22.70 1368.00 1.0000000 -10530 22.80 1368.00 1.0000000 -10520 29.70 1368.00 1.0000000 -10510 44.80 1368.00 1.0000000 -10500 57.70 1368.00 1.0000000 -10490 48.50 1368.00 1.0000000 -10480 34.60 1368.00 1.0000000 -10470 32.80 1368.00 1.0000000 -10460 30.90 1368.00 1.0000000 -10450 28.80 1368.00 1.0000000 -10440 31.80 1368.00 1.0000000 -10430 33.00 1368.00 1.0000000 -10420 29.60 1368.00 1.0000000 -10410 35.30 1368.00 1.0000000 -10400 49.90 1368.00 1.0000000 -10390 52.70 1368.00 1.0000000 -10380 52.80 1368.00 1.0000000 -10370 59.00 1368.00 1.0000000 -10360 64.30 1368.00 1.0000000 -10350 61.90 1368.00 1.0000000 -10340 50.50 1368.00 1.0000000 -10330 47.20 1368.00 1.0000000 -10320 49.50 1368.00 1.0000000 -10310 42.60 1368.00 1.0000000 -10300 26.10 1368.00 1.0000000 -10290 15.60 1368.00 1.0000000 -10280 11.50 1368.00 1.0000000 -10270 8.70 1368.00 1.0000000 -10260 10.30 1368.00 1.0000000 -10250 10.40 1368.00 1.0000000 -10240 7.50 1368.00 1.0000000 -10230 8.90 1368.00 1.0000000 -10220 9.00 1368.00 1.0000000 -10210 4.80 1368.00 1.0000000 -10200 5.70 1368.00 1.0000000 -10190 13.70 1368.00 1.0000000 -10180 16.00 1368.00 1.0000000 -10170 8.70 1368.00 1.0000000 -10160 6.60 1368.00 1.0000000 -10150 17.40 1368.00 1.0000000 -10140 32.00 1368.00 1.0000000 -10130 36.00 1368.00 1.0000000 -10120 37.10 1368.00 1.0000000 -10110 44.40 1368.00 1.0000000 -10100 42.50 1368.00 1.0000000 -10090 33.20 1368.00 1.0000000 -10080 33.80 1368.00 1.0000000 -10070 37.00 1368.00 1.0000000 -10060 42.90 1368.00 1.0000000 -10050 61.30 1368.00 1.0000000 -10040 67.80 1368.00 1.0000000 -10030 53.70 1368.00 1.0000000 -10020 45.90 1368.00 1.0000000 -10010 46.00 1368.00 1.0000000 -10000 50.60 1368.00 1.0000000 -9990 51.30 1368.00 1.0000000 -9980 37.00 1368.00 1.0000000 -9970 23.60 1368.00 1.0000000 -9960 22.30 1368.00 1.0000000 -9950 23.00 1368.00 1.0000000 -9940 26.20 1368.00 1.0000000 -9930 28.80 1368.00 1.0000000 -9920 27.00 1368.00 1.0000000 -9910 36.80 1368.00 1.0000000 -9900 50.10 1368.00 1.0000000 -9890 46.30 1368.00 1.0000000 -9880 39.70 1368.00 1.0000000 -9870 47.40 1368.00 1.0000000 -9860 58.90 1368.00 1.0000000 -9850 53.90 1368.00 1.0000000 -9840 44.50 1368.00 1.0000000 -9830 42.00 1368.00 1.0000000 -9820 44.40 1368.00 1.0000000 -9810 47.30 1368.00 1.0000000 -9800 49.00 1368.00 1.0000000 -9790 50.90 1368.00 1.0000000 -9780 52.10 1368.00 1.0000000 -9770 48.20 1368.00 1.0000000 -9760 39.90 1368.00 1.0000000 -9750 34.40 1368.00 1.0000000 -9740 33.80 1368.00 1.0000000 -9730 35.90 1368.00 1.0000000 -9720 39.80 1368.00 1.0000000 -9710 48.70 1368.00 1.0000000 -9700 57.30 1368.00 1.0000000 -9690 51.70 1368.00 1.0000000 -9680 48.10 1368.00 1.0000000 -9670 51.90 1368.00 1.0000000 -9660 48.70 1368.00 1.0000000 -9650 48.90 1368.00 1.0000000 -9640 57.20 1368.00 1.0000000 -9630 61.00 1368.00 1.0000000 -9620 49.30 1368.00 1.0000000 -9610 27.80 1368.00 1.0000000 -9600 13.20 1368.00 1.0000000 -9590 5.50 1368.00 1.0000000 -9580 4.50 1368.00 1.0000000 -9570 12.80 1368.00 1.0000000 -9560 17.70 1368.00 1.0000000 -9550 12.80 1368.00 1.0000000 -9540 8.10 1368.00 1.0000000 -9530 3.10 1368.00 1.0000000 -9520 3.60 1368.00 1.0000000 -9510 16.30 1368.00 1.0000000 -9500 23.20 1368.00 1.0000000 -9490 13.30 1368.00 1.0000000 -9480 8.20 1368.00 1.0000000 -9470 10.50 1368.00 1.0000000 -9460 7.50 1368.00 1.0000000 -9450 7.40 1368.00 1.0000000 -9440 14.90 1368.00 1.0000000 -9430 18.30 1368.00 1.0000000 -9420 22.00 1368.00 1.0000000 -9410 32.30 1368.00 1.0000000 -9400 44.20 1368.00 1.0000000 -9390 44.20 1368.00 1.0000000 -9380 33.70 1368.00 1.0000000 -9370 31.60 1368.00 1.0000000 -9360 29.30 1368.00 1.0000000 -9350 25.20 1368.00 1.0000000 -9340 23.70 1368.00 1.0000000 -9330 20.20 1368.00 1.0000000 -9320 16.50 1368.00 1.0000000 -9310 14.00 1368.00 1.0000000 -9300 12.00 1368.00 1.0000000 -9290 14.90 1368.00 1.0000000 -9280 26.60 1368.00 1.0000000 -9270 36.70 1368.00 1.0000000 -9260 35.00 1368.00 1.0000000 -9250 30.80 1368.00 1.0000000 -9240 38.00 1368.00 1.0000000 -9230 63.50 1368.00 1.0000000 -9220 76.10 1368.00 1.0000000 -9210 56.40 1368.00 1.0000000 -9200 39.00 1368.00 1.0000000 -9190 25.10 1368.00 1.0000000 -9180 21.10 1368.00 1.0000000 -9170 31.70 1368.00 1.0000000 -9160 42.70 1368.00 1.0000000 -9150 37.90 1368.00 1.0000000 -9140 31.30 1368.00 1.0000000 -9130 33.50 1368.00 1.0000000 -9120 38.80 1368.00 1.0000000 -9110 50.00 1368.00 1.0000000 -9100 53.60 1368.00 1.0000000 -9090 34.60 1368.00 1.0000000 -9080 19.50 1368.00 1.0000000 -9070 12.60 1368.00 1.0000000 -9060 7.50 1368.00 1.0000000 -9050 5.40 1368.00 1.0000000 -9040 4.90 1368.00 1.0000000 -9030 7.60 1368.00 1.0000000 -9020 13.50 1368.00 1.0000000 -9010 17.60 1368.00 1.0000000 -9000 19.20 1368.00 1.0000000 -8990 20.30 1368.00 1.0000000 -8980 29.00 1368.00 1.0000000 -8970 39.90 1368.00 1.0000000 -8960 34.30 1368.00 1.0000000 -8950 31.60 1368.00 1.0000000 -8940 44.70 1368.00 1.0000000 -8930 47.50 1368.00 1.0000000 -8920 36.00 1368.00 1.0000000 -8910 30.90 1368.00 1.0000000 -8900 38.00 1368.00 1.0000000 -8890 63.20 1368.00 1.0000000 -8880 81.80 1368.00 1.0000000 -8870 66.30 1368.00 1.0000000 -8860 59.30 1368.00 1.0000000 -8850 58.00 1368.00 1.0000000 -8840 46.20 1368.00 1.0000000 -8830 37.20 1368.00 1.0000000 -8820 38.20 1368.00 1.0000000 -8810 42.50 1368.00 1.0000000 -8800 49.20 1368.00 1.0000000 -8790 55.50 1368.00 1.0000000 -8780 53.70 1368.00 1.0000000 -8770 50.30 1368.00 1.0000000 -8760 48.20 1368.00 1.0000000 -8750 47.20 1368.00 1.0000000 -8740 53.10 1368.00 1.0000000 -8730 65.40 1368.00 1.0000000 -8720 69.80 1368.00 1.0000000 -8710 57.40 1368.00 1.0000000 -8700 44.00 1368.00 1.0000000 -8690 40.00 1368.00 1.0000000 -8680 45.90 1368.00 1.0000000 -8670 48.20 1368.00 1.0000000 -8660 35.30 1368.00 1.0000000 -8650 22.40 1368.00 1.0000000 -8640 23.10 1368.00 1.0000000 -8630 34.30 1368.00 1.0000000 -8620 42.00 1368.00 1.0000000 -8610 38.00 1368.00 1.0000000 -8600 30.40 1368.00 1.0000000 -8590 30.60 1368.00 1.0000000 -8580 34.40 1368.00 1.0000000 -8570 35.70 1368.00 1.0000000 -8560 44.20 1368.00 1.0000000 -8550 53.70 1368.00 1.0000000 -8540 51.10 1368.00 1.0000000 -8530 49.30 1368.00 1.0000000 -8520 48.10 1368.00 1.0000000 -8510 41.00 1368.00 1.0000000 -8500 40.10 1368.00 1.0000000 -8490 41.50 1368.00 1.0000000 -8480 33.50 1368.00 1.0000000 -8470 23.90 1368.00 1.0000000 -8460 19.60 1368.00 1.0000000 -8450 17.30 1368.00 1.0000000 -8440 14.90 1368.00 1.0000000 -8430 14.60 1368.00 1.0000000 -8420 15.30 1368.00 1.0000000 -8410 15.00 1368.00 1.0000000 -8400 13.00 1368.00 1.0000000 -8390 10.80 1368.00 1.0000000 -8380 12.20 1368.00 1.0000000 -8370 18.00 1368.00 1.0000000 -8360 22.50 1368.00 1.0000000 -8350 21.90 1368.00 1.0000000 -8340 20.80 1368.00 1.0000000 -8330 27.70 1368.00 1.0000000 -8320 41.10 1368.00 1.0000000 -8310 46.00 1368.00 1.0000000 -8300 47.30 1368.00 1.0000000 -8290 54.30 1368.00 1.0000000 -8280 51.30 1368.00 1.0000000 -8270 40.90 1368.00 1.0000000 -8260 34.70 1368.00 1.0000000 -8250 28.60 1368.00 1.0000000 -8240 19.60 1368.00 1.0000000 -8230 12.20 1368.00 1.0000000 -8220 9.60 1368.00 1.0000000 -8210 11.80 1368.00 1.0000000 -8200 19.70 1368.00 1.0000000 -8190 28.50 1368.00 1.0000000 -8180 28.40 1368.00 1.0000000 -8170 18.20 1368.00 1.0000000 -8160 15.40 1368.00 1.0000000 -8150 34.90 1368.00 1.0000000 -8140 61.20 1368.00 1.0000000 -8130 61.70 1368.00 1.0000000 -8120 47.90 1368.00 1.0000000 -8110 39.30 1368.00 1.0000000 -8100 33.10 1368.00 1.0000000 -8090 25.80 1368.00 1.0000000 -8080 21.80 1368.00 1.0000000 -8070 20.10 1368.00 1.0000000 -8060 22.60 1368.00 1.0000000 -8050 31.70 1368.00 1.0000000 -8040 32.40 1368.00 1.0000000 -8030 22.20 1368.00 1.0000000 -8020 17.70 1368.00 1.0000000 -8010 17.20 1368.00 1.0000000 -8000 12.50 1368.00 1.0000000 -7990 6.10 1368.00 1.0000000 -7980 5.10 1368.00 1.0000000 -7970 11.30 1368.00 1.0000000 -7960 28.10 1368.00 1.0000000 -7950 48.40 1368.00 1.0000000 -7940 44.60 1368.00 1.0000000 -7930 26.50 1368.00 1.0000000 -7920 21.00 1368.00 1.0000000 -7910 22.90 1368.00 1.0000000 -7900 22.30 1368.00 1.0000000 -7890 20.50 1368.00 1.0000000 -7880 22.00 1368.00 1.0000000 -7870 28.20 1368.00 1.0000000 -7860 30.80 1368.00 1.0000000 -7850 21.20 1368.00 1.0000000 -7840 17.90 1368.00 1.0000000 -7830 27.00 1368.00 1.0000000 -7820 28.70 1368.00 1.0000000 -7810 21.70 1368.00 1.0000000 -7800 20.50 1368.00 1.0000000 -7790 22.70 1368.00 1.0000000 -7780 25.10 1368.00 1.0000000 -7770 28.80 1368.00 1.0000000 -7760 30.90 1368.00 1.0000000 -7750 26.40 1368.00 1.0000000 -7740 17.60 1368.00 1.0000000 -7730 9.70 1368.00 1.0000000 -7720 6.70 1368.00 1.0000000 -7710 9.50 1368.00 1.0000000 -7700 16.70 1368.00 1.0000000 -7690 23.50 1368.00 1.0000000 -7680 21.80 1368.00 1.0000000 -7670 15.70 1368.00 1.0000000 -7660 16.40 1368.00 1.0000000 -7650 17.10 1368.00 1.0000000 -7640 10.50 1368.00 1.0000000 -7630 2.70 1368.00 1.0000000 -7620 -2.30 1368.00 1.0000000 -7610 0.10 1368.00 1.0000000 -7600 12.60 1368.00 1.0000000 -7590 28.30 1368.00 1.0000000 -7580 29.90 1368.00 1.0000000 -7570 17.20 1368.00 1.0000000 -7560 11.60 1368.00 1.0000000 -7550 14.60 1368.00 1.0000000 -7540 15.70 1368.00 1.0000000 -7530 14.50 1368.00 1.0000000 -7520 14.90 1368.00 1.0000000 -7510 19.20 1368.00 1.0000000 -7500 19.60 1368.00 1.0000000 -7490 6.20 1368.00 1.0000000 -7480 -6.40 1368.00 1.0000000 -7470 -7.40 1368.00 1.0000000 -7460 0.70 1368.00 1.0000000 -7450 12.20 1368.00 1.0000000 -7440 16.30 1368.00 1.0000000 -7430 11.30 1368.00 1.0000000 -7420 11.20 1368.00 1.0000000 -7410 21.40 1368.00 1.0000000 -7400 26.10 1368.00 1.0000000 -7390 18.70 1368.00 1.0000000 -7380 11.80 1368.00 1.0000000 -7370 9.50 1368.00 1.0000000 -7360 12.20 1368.00 1.0000000 -7350 15.90 1368.00 1.0000000 -7340 16.70 1368.00 1.0000000 -7330 13.00 1368.00 1.0000000 -7320 7.40 1368.00 1.0000000 -7310 3.00 1368.00 1.0000000 -7300 1.50 1368.00 1.0000000 -7290 4.50 1368.00 1.0000000 -7280 8.50 1368.00 1.0000000 -7270 7.80 1368.00 1.0000000 -7260 7.10 1368.00 1.0000000 -7250 10.80 1368.00 1.0000000 -7240 13.40 1368.00 1.0000000 -7230 13.50 1368.00 1.0000000 -7220 8.20 1368.00 1.0000000 -7210 1.00 1368.00 1.0000000 -7200 2.90 1368.00 1.0000000 -7190 13.90 1368.00 1.0000000 -7180 28.30 1368.00 1.0000000 -7170 37.40 1368.00 1.0000000 -7160 30.50 1368.00 1.0000000 -7150 18.20 1368.00 1.0000000 -7140 18.10 1368.00 1.0000000 -7130 27.10 1368.00 1.0000000 -7120 34.20 1368.00 1.0000000 -7110 36.40 1368.00 1.0000000 -7100 34.60 1368.00 1.0000000 -7090 28.60 1368.00 1.0000000 -7080 21.80 1368.00 1.0000000 -7070 18.20 1368.00 1.0000000 -7060 16.10 1368.00 1.0000000 -7050 15.20 1368.00 1.0000000 -7040 16.70 1368.00 1.0000000 -7030 20.40 1368.00 1.0000000 -7020 24.30 1368.00 1.0000000 -7010 23.40 1368.00 1.0000000 -7000 20.70 1368.00 1.0000000 -6990 22.70 1368.00 1.0000000 -6980 26.90 1368.00 1.0000000 -6970 27.20 1368.00 1.0000000 -6960 23.60 1368.00 1.0000000 -6950 19.50 1368.00 1.0000000 -6940 18.90 1368.00 1.0000000 -6930 20.60 1368.00 1.0000000 -6920 18.50 1368.00 1.0000000 -6910 15.70 1368.00 1.0000000 -6900 20.60 1368.00 1.0000000 -6890 32.60 1368.00 1.0000000 -6880 45.00 1368.00 1.0000000 -6870 41.20 1368.00 1.0000000 -6860 23.60 1368.00 1.0000000 -6850 15.30 1368.00 1.0000000 -6840 19.00 1368.00 1.0000000 -6830 29.40 1368.00 1.0000000 -6820 30.80 1368.00 1.0000000 -6810 21.30 1368.00 1.0000000 -6800 17.90 1368.00 1.0000000 -6790 15.10 1368.00 1.0000000 -6780 9.40 1368.00 1.0000000 -6770 12.90 1368.00 1.0000000 -6760 28.10 1368.00 1.0000000 -6750 37.10 1368.00 1.0000000 -6740 29.10 1368.00 1.0000000 -6730 18.60 1368.00 1.0000000 -6720 13.60 1368.00 1.0000000 -6710 14.50 1368.00 1.0000000 -6700 14.60 1368.00 1.0000000 -6690 12.40 1368.00 1.0000000 -6680 13.50 1368.00 1.0000000 -6670 17.60 1368.00 1.0000000 -6660 22.00 1368.00 1.0000000 -6650 27.90 1368.00 1.0000000 -6640 33.60 1368.00 1.0000000 -6630 30.40 1368.00 1.0000000 -6620 24.90 1368.00 1.0000000 -6610 24.50 1368.00 1.0000000 -6600 25.20 1368.00 1.0000000 -6590 26.60 1368.00 1.0000000 -6580 30.90 1368.00 1.0000000 -6570 33.50 1368.00 1.0000000 -6560 27.30 1368.00 1.0000000 -6550 20.20 1368.00 1.0000000 -6540 20.20 1368.00 1.0000000 -6530 23.60 1368.00 1.0000000 -6520 26.60 1368.00 1.0000000 -6510 24.30 1368.00 1.0000000 -6500 19.60 1368.00 1.0000000 -6490 23.30 1368.00 1.0000000 -6480 29.10 1368.00 1.0000000 -6470 21.00 1368.00 1.0000000 -6460 9.60 1368.00 1.0000000 -6450 8.10 1368.00 1.0000000 -6440 14.70 1368.00 1.0000000 -6430 27.70 1368.00 1.0000000 -6420 35.40 1368.00 1.0000000 -6410 29.00 1368.00 1.0000000 -6400 26.50 1368.00 1.0000000 -6390 32.70 1368.00 1.0000000 -6380 32.60 1368.00 1.0000000 -6370 22.40 1368.00 1.0000000 -6360 13.50 1368.00 1.0000000 -6350 7.50 1368.00 1.0000000 -6340 2.00 1368.00 1.0000000 -6330 -0.50 1368.00 1.0000000 -6320 2.60 1368.00 1.0000000 -6310 10.70 1368.00 1.0000000 -6300 20.60 1368.00 1.0000000 -6290 29.50 1368.00 1.0000000 -6280 34.00 1368.00 1.0000000 -6270 27.70 1368.00 1.0000000 -6260 16.30 1368.00 1.0000000 -6250 12.50 1368.00 1.0000000 -6240 10.50 1368.00 1.0000000 -6230 2.10 1368.00 1.0000000 -6220 0.10 1368.00 1.0000000 -6210 13.60 1368.00 1.0000000 -6200 31.60 1368.00 1.0000000 -6190 39.70 1368.00 1.0000000 -6180 31.90 1368.00 1.0000000 -6170 19.80 1368.00 1.0000000 -6160 17.50 1368.00 1.0000000 -6150 20.40 1368.00 1.0000000 -6140 27.00 1368.00 1.0000000 -6130 38.70 1368.00 1.0000000 -6120 40.10 1368.00 1.0000000 -6110 36.90 1368.00 1.0000000 -6100 41.60 1368.00 1.0000000 -6090 43.10 1368.00 1.0000000 -6080 42.30 1368.00 1.0000000 -6070 43.10 1368.00 1.0000000 -6060 31.50 1368.00 1.0000000 -6050 14.60 1368.00 1.0000000 -6040 13.10 1368.00 1.0000000 -6030 26.90 1368.00 1.0000000 -6020 37.80 1368.00 1.0000000 -6010 32.10 1368.00 1.0000000 -6000 23.20 1368.00 1.0000000 -5990 17.30 1368.00 1.0000000 -5980 12.10 1368.00 1.0000000 -5970 7.90 1368.00 1.0000000 -5960 5.80 1368.00 1.0000000 -5950 4.90 1368.00 1.0000000 -5940 5.90 1368.00 1.0000000 -5930 9.40 1368.00 1.0000000 -5920 11.60 1368.00 1.0000000 -5910 12.10 1368.00 1.0000000 -5900 15.00 1368.00 1.0000000 -5890 23.80 1368.00 1.0000000 -5880 33.80 1368.00 1.0000000 -5870 38.30 1368.00 1.0000000 -5860 41.10 1368.00 1.0000000 -5850 44.00 1368.00 1.0000000 -5840 44.70 1368.00 1.0000000 -5830 41.10 1368.00 1.0000000 -5820 33.40 1368.00 1.0000000 -5810 23.90 1368.00 1.0000000 -5800 17.80 1368.00 1.0000000 -5790 17.10 1368.00 1.0000000 -5780 17.60 1368.00 1.0000000 -5770 18.10 1368.00 1.0000000 -5760 20.90 1368.00 1.0000000 -5750 25.90 1368.00 1.0000000 -5740 30.60 1368.00 1.0000000 -5730 31.60 1368.00 1.0000000 -5720 24.60 1368.00 1.0000000 -5710 16.70 1368.00 1.0000000 -5700 17.30 1368.00 1.0000000 -5690 23.80 1368.00 1.0000000 -5680 26.90 1368.00 1.0000000 -5670 21.40 1368.00 1.0000000 -5660 12.00 1368.00 1.0000000 -5650 2.30 1368.00 1.0000000 -5640 -4.40 1368.00 1.0000000 -5630 -4.70 1368.00 1.0000000 -5620 0.90 1368.00 1.0000000 -5610 9.10 1368.00 1.0000000 -5600 16.80 1368.00 1.0000000 -5590 23.60 1368.00 1.0000000 -5580 28.00 1368.00 1.0000000 -5570 33.40 1368.00 1.0000000 -5560 38.30 1368.00 1.0000000 -5550 31.10 1368.00 1.0000000 -5540 19.60 1368.00 1.0000000 -5530 12.90 1368.00 1.0000000 -5520 9.00 1368.00 1.0000000 -5510 5.50 1368.00 1.0000000 -5500 3.10 1368.00 1.0000000 -5490 6.50 1368.00 1.0000000 -5480 16.40 1368.00 1.0000000 -5470 26.80 1368.00 1.0000000 -5460 31.10 1368.00 1.0000000 -5450 33.10 1368.00 1.0000000 -5440 39.90 1368.00 1.0000000 -5430 47.10 1368.00 1.0000000 -5420 50.40 1368.00 1.0000000 -5410 53.20 1368.00 1.0000000 -5400 48.30 1368.00 1.0000000 -5390 30.10 1368.00 1.0000000 -5380 13.20 1368.00 1.0000000 -5370 5.60 1368.00 1.0000000 -5360 3.20 1368.00 1.0000000 -5350 2.30 1368.00 1.0000000 -5340 2.70 1368.00 1.0000000 -5330 4.80 1368.00 1.0000000 -5320 7.70 1368.00 1.0000000 -5310 10.50 1368.00 1.0000000 -5300 13.90 1368.00 1.0000000 -5290 17.10 1368.00 1.0000000 -5280 18.30 1368.00 1.0000000 -5270 17.30 1368.00 1.0000000 -5260 19.20 1368.00 1.0000000 -5250 28.00 1368.00 1.0000000 -5240 37.60 1368.00 1.0000000 -5230 44.50 1368.00 1.0000000 -5220 52.60 1368.00 1.0000000 -5210 53.40 1368.00 1.0000000 -5200 43.80 1368.00 1.0000000 -5190 35.40 1368.00 1.0000000 -5180 33.80 1368.00 1.0000000 -5170 44.80 1368.00 1.0000000 -5160 62.20 1368.00 1.0000000 -5150 61.80 1368.00 1.0000000 -5140 50.20 1368.00 1.0000000 -5130 43.50 1368.00 1.0000000 -5120 38.80 1368.00 1.0000000 -5110 30.10 1368.00 1.0000000 -5100 21.00 1368.00 1.0000000 -5090 19.30 1368.00 1.0000000 -5080 27.20 1368.00 1.0000000 -5070 40.40 1368.00 1.0000000 -5060 46.40 1368.00 1.0000000 -5050 42.20 1368.00 1.0000000 -5040 34.10 1368.00 1.0000000 -5030 26.10 1368.00 1.0000000 -5020 23.30 1368.00 1.0000000 -5010 26.60 1368.00 1.0000000 -5000 34.20 1368.00 1.0000000 -4990 43.10 1368.00 1.0000000 -4980 46.70 1368.00 1.0000000 -4970 46.20 1368.00 1.0000000 -4960 49.90 1368.00 1.0000000 -4950 53.40 1368.00 1.0000000 -4940 46.10 1368.00 1.0000000 -4930 31.40 1368.00 1.0000000 -4920 21.40 1368.00 1.0000000 -4910 18.20 1368.00 1.0000000 -4900 14.00 1368.00 1.0000000 -4890 6.40 1368.00 1.0000000 -4880 1.30 1368.00 1.0000000 -4870 0.50 1368.00 1.0000000 -4860 2.30 1368.00 1.0000000 -4850 4.40 1368.00 1.0000000 -4840 6.80 1368.00 1.0000000 -4830 16.60 1368.00 1.0000000 -4820 36.60 1368.00 1.0000000 -4810 51.80 1368.00 1.0000000 -4800 47.80 1368.00 1.0000000 -4790 39.20 1368.00 1.0000000 -4780 37.80 1368.00 1.0000000 -4770 36.20 1368.00 1.0000000 -4760 31.00 1368.00 1.0000000 -4750 30.50 1368.00 1.0000000 -4740 37.70 1368.00 1.0000000 -4730 51.40 1368.00 1.0000000 -4720 54.50 1368.00 1.0000000 -4710 41.60 1368.00 1.0000000 -4700 36.20 1368.00 1.0000000 -4690 38.40 1368.00 1.0000000 -4680 36.50 1368.00 1.0000000 -4670 33.40 1368.00 1.0000000 -4660 39.70 1368.00 1.0000000 -4650 48.30 1368.00 1.0000000 -4640 42.70 1368.00 1.0000000 -4630 32.50 1368.00 1.0000000 -4620 34.70 1368.00 1.0000000 -4610 42.10 1368.00 1.0000000 -4600 40.40 1368.00 1.0000000 -4590 31.60 1368.00 1.0000000 -4580 22.70 1368.00 1.0000000 -4570 19.40 1368.00 1.0000000 -4560 21.90 1368.00 1.0000000 -4550 29.40 1368.00 1.0000000 -4540 40.50 1368.00 1.0000000 -4530 53.90 1368.00 1.0000000 -4520 61.60 1368.00 1.0000000 -4510 48.40 1368.00 1.0000000 -4500 34.20 1368.00 1.0000000 -4490 31.20 1368.00 1.0000000 -4480 24.10 1368.00 1.0000000 -4470 14.10 1368.00 1.0000000 -4460 14.00 1368.00 1.0000000 -4450 20.30 1368.00 1.0000000 -4440 26.00 1368.00 1.0000000 -4430 31.70 1368.00 1.0000000 -4420 33.10 1368.00 1.0000000 -4410 27.50 1368.00 1.0000000 -4400 27.70 1368.00 1.0000000 -4390 38.80 1368.00 1.0000000 -4380 49.40 1368.00 1.0000000 -4370 49.50 1368.00 1.0000000 -4360 42.10 1368.00 1.0000000 -4350 35.90 1368.00 1.0000000 -4340 37.10 1368.00 1.0000000 -4330 42.50 1368.00 1.0000000 -4320 43.90 1368.00 1.0000000 -4310 38.90 1368.00 1.0000000 -4300 31.90 1368.00 1.0000000 -4290 27.00 1368.00 1.0000000 -4280 27.00 1368.00 1.0000000 -4270 35.20 1368.00 1.0000000 -4260 49.70 1368.00 1.0000000 -4250 55.10 1368.00 1.0000000 -4240 48.70 1368.00 1.0000000 -4230 49.70 1368.00 1.0000000 -4220 47.40 1368.00 1.0000000 -4210 31.90 1368.00 1.0000000 -4200 24.50 1368.00 1.0000000 -4190 31.40 1368.00 1.0000000 -4180 44.20 1368.00 1.0000000 -4170 54.00 1368.00 1.0000000 -4160 49.10 1368.00 1.0000000 -4150 32.50 1368.00 1.0000000 -4140 22.90 1368.00 1.0000000 -4130 23.20 1368.00 1.0000000 -4120 27.30 1368.00 1.0000000 -4110 35.10 1368.00 1.0000000 -4100 46.90 1368.00 1.0000000 -4090 59.40 1368.00 1.0000000 -4080 64.90 1368.00 1.0000000 -4070 60.40 1368.00 1.0000000 -4060 54.30 1368.00 1.0000000 -4050 43.70 1368.00 1.0000000 -4040 32.30 1368.00 1.0000000 -4030 28.00 1368.00 1.0000000 -4020 30.40 1368.00 1.0000000 -4010 41.80 1368.00 1.0000000 -4000 54.10 1368.00 1.0000000 -3990 47.90 1368.00 1.0000000 -3980 39.40 1368.00 1.0000000 -3970 40.60 1368.00 1.0000000 -3960 39.70 1368.00 1.0000000 -3950 36.40 1368.00 1.0000000 -3940 36.50 1368.00 1.0000000 -3930 36.90 1368.00 1.0000000 -3920 37.30 1368.00 1.0000000 -3910 40.40 1368.00 1.0000000 -3900 36.60 1368.00 1.0000000 -3890 24.10 1368.00 1.0000000 -3880 17.50 1368.00 1.0000000 -3870 19.40 1368.00 1.0000000 -3860 29.40 1368.00 1.0000000 -3850 48.20 1368.00 1.0000000 -3840 54.70 1368.00 1.0000000 -3830 43.40 1368.00 1.0000000 -3820 40.20 1368.00 1.0000000 -3810 49.10 1368.00 1.0000000 -3800 58.80 1368.00 1.0000000 -3790 51.00 1368.00 1.0000000 -3780 38.80 1368.00 1.0000000 -3770 38.10 1368.00 1.0000000 -3760 36.60 1368.00 1.0000000 -3750 26.30 1368.00 1.0000000 -3740 21.30 1368.00 1.0000000 -3730 29.50 1368.00 1.0000000 -3720 43.20 1368.00 1.0000000 -3710 45.10 1368.00 1.0000000 -3700 29.70 1368.00 1.0000000 -3690 18.40 1368.00 1.0000000 -3680 22.70 1368.00 1.0000000 -3670 31.70 1368.00 1.0000000 -3660 36.20 1368.00 1.0000000 -3650 37.20 1368.00 1.0000000 -3640 33.20 1368.00 1.0000000 -3630 28.50 1368.00 1.0000000 -3620 25.50 1368.00 1.0000000 -3610 23.80 1368.00 1.0000000 -3600 26.40 1368.00 1.0000000 -3590 32.70 1368.00 1.0000000 -3580 39.30 1368.00 1.0000000 -3570 46.10 1368.00 1.0000000 -3560 50.60 1368.00 1.0000000 -3550 48.40 1368.00 1.0000000 -3540 36.60 1368.00 1.0000000 -3530 21.40 1368.00 1.0000000 -3520 16.10 1368.00 1.0000000 -3510 18.10 1368.00 1.0000000 -3500 19.70 1368.00 1.0000000 -3490 26.10 1368.00 1.0000000 -3480 35.50 1368.00 1.0000000 -3470 36.10 1368.00 1.0000000 -3460 31.10 1368.00 1.0000000 -3450 26.60 1368.00 1.0000000 -3440 22.50 1368.00 1.0000000 -3430 20.50 1368.00 1.0000000 -3420 18.20 1368.00 1.0000000 -3410 15.20 1368.00 1.0000000 -3400 14.40 1368.00 1.0000000 -3390 13.20 1368.00 1.0000000 -3380 12.20 1368.00 1.0000000 -3370 15.90 1368.00 1.0000000 -3360 23.50 1368.00 1.0000000 -3350 34.30 1368.00 1.0000000 -3340 42.80 1368.00 1.0000000 -3330 34.00 1368.00 1.0000000 -3320 23.00 1368.00 1.0000000 -3310 22.60 1368.00 1.0000000 -3300 24.60 1368.00 1.0000000 -3290 29.90 1368.00 1.0000000 -3280 32.80 1368.00 1.0000000 -3270 22.90 1368.00 1.0000000 -3260 17.40 1368.00 1.0000000 -3250 26.90 1368.00 1.0000000 -3240 38.80 1368.00 1.0000000 -3230 35.10 1368.00 1.0000000 -3220 25.40 1368.00 1.0000000 -3210 25.40 1368.00 1.0000000 -3200 27.70 1368.00 1.0000000 -3190 31.90 1368.00 1.0000000 -3180 39.80 1368.00 1.0000000 -3170 38.30 1368.00 1.0000000 -3160 37.50 1368.00 1.0000000 -3150 44.60 1368.00 1.0000000 -3140 38.10 1368.00 1.0000000 -3130 26.40 1368.00 1.0000000 -3120 27.40 1368.00 1.0000000 -3110 34.30 1368.00 1.0000000 -3100 39.20 1368.00 1.0000000 -3090 40.50 1368.00 1.0000000 -3080 42.30 1368.00 1.0000000 -3070 44.80 1368.00 1.0000000 -3060 41.00 1368.00 1.0000000 -3050 36.80 1368.00 1.0000000 -3040 39.10 1368.00 1.0000000 -3030 40.70 1368.00 1.0000000 -3020 35.60 1368.00 1.0000000 -3010 29.00 1368.00 1.0000000 -3000 28.70 1368.00 1.0000000 -2990 30.90 1368.00 1.0000000 -2980 29.30 1368.00 1.0000000 -2970 35.00 1368.00 1.0000000 -2960 46.10 1368.00 1.0000000 -2950 44.50 1368.00 1.0000000 -2940 36.70 1368.00 1.0000000 -2930 32.10 1368.00 1.0000000 -2920 28.50 1368.00 1.0000000 -2910 23.70 1368.00 1.0000000 -2900 21.50 1368.00 1.0000000 -2890 28.50 1368.00 1.0000000 -2880 38.40 1368.00 1.0000000 -2870 41.90 1368.00 1.0000000 -2860 42.30 1368.00 1.0000000 -2850 36.60 1368.00 1.0000000 -2840 25.10 1368.00 1.0000000 -2830 19.00 1368.00 1.0000000 -2820 16.20 1368.00 1.0000000 -2810 8.30 1368.00 1.0000000 -2800 1.60 1368.00 1.0000000 -2790 2.30 1368.00 1.0000000 -2780 3.80 1368.00 1.0000000 -2770 -0.80 1368.00 1.0000000 -2760 -3.30 1368.00 1.0000000 -2750 1.50 1368.00 1.0000000 -2740 7.50 1368.00 1.0000000 -2730 12.00 1368.00 1.0000000 -2720 16.10 1368.00 1.0000000 -2710 20.80 1368.00 1.0000000 -2700 29.70 1368.00 1.0000000 -2690 40.80 1368.00 1.0000000 -2680 39.20 1368.00 1.0000000 -2670 28.80 1368.00 1.0000000 -2660 29.40 1368.00 1.0000000 -2650 38.20 1368.00 1.0000000 -2640 44.80 1368.00 1.0000000 -2630 47.80 1368.00 1.0000000 -2620 51.40 1368.00 1.0000000 -2610 49.70 1368.00 1.0000000 -2600 42.80 1368.00 1.0000000 -2590 44.80 1368.00 1.0000000 -2580 52.20 1368.00 1.0000000 -2570 51.20 1368.00 1.0000000 -2560 43.40 1368.00 1.0000000 -2550 35.60 1368.00 1.0000000 -2540 33.60 1368.00 1.0000000 -2530 33.70 1368.00 1.0000000 -2520 32.40 1368.00 1.0000000 -2510 35.00 1368.00 1.0000000 -2500 38.20 1368.00 1.0000000 -2490 38.70 1368.00 1.0000000 -2480 44.30 1368.00 1.0000000 -2470 53.30 1368.00 1.0000000 -2460 54.50 1368.00 1.0000000 -2450 55.70 1368.00 1.0000000 -2440 60.60 1368.00 1.0000000 -2430 56.20 1368.00 1.0000000 -2420 36.90 1368.00 1.0000000 -2410 16.10 1368.00 1.0000000 -2400 6.60 1368.00 1.0000000 -2390 6.10 1368.00 1.0000000 -2380 8.40 1368.00 1.0000000 -2370 7.80 1368.00 1.0000000 -2360 3.90 1368.00 1.0000000 -2350 4.00 1368.00 1.0000000 -2340 10.70 1368.00 1.0000000 -2330 20.90 1368.00 1.0000000 -2320 28.80 1368.00 1.0000000 -2310 33.50 1368.00 1.0000000 -2300 41.00 1368.00 1.0000000 -2290 50.80 1368.00 1.0000000 -2280 54.90 1368.00 1.0000000 -2270 47.40 1368.00 1.0000000 -2260 42.10 1368.00 1.0000000 -2250 42.90 1368.00 1.0000000 -2240 38.40 1368.00 1.0000000 -2230 38.60 1368.00 1.0000000 -2220 37.30 1368.00 1.0000000 -2210 25.70 1368.00 1.0000000 -2200 23.50 1368.00 1.0000000 -2190 26.90 1368.00 1.0000000 -2180 22.80 1368.00 1.0000000 -2170 19.60 1368.00 1.0000000 -2160 22.40 1368.00 1.0000000 -2150 28.40 1368.00 1.0000000 -2140 35.50 1368.00 1.0000000 -2130 36.80 1368.00 1.0000000 -2120 33.40 1368.00 1.0000000 -2110 31.90 1368.00 1.0000000 -2100 30.90 1368.00 1.0000000 -2090 34.00 1368.00 1.0000000 -2080 40.70 1368.00 1.0000000 -2070 39.80 1368.00 1.0000000 -2060 30.10 1368.00 1.0000000 -2050 23.00 1368.00 1.0000000 -2040 24.30 1368.00 1.0000000 -2030 30.20 1368.00 1.0000000 -2020 34.60 1368.00 1.0000000 -2010 32.60 1368.00 1.0000000 -2000 28.60 1368.00 1.0000000 -1990 26.80 1368.00 1.0000000 -1980 25.10 1368.00 1.0000000 -1970 29.20 1368.00 1.0000000 -1960 32.50 1368.00 1.0000000 -1950 28.30 1368.00 1.0000000 -1940 24.00 1368.00 1.0000000 -1930 19.80 1368.00 1.0000000 -1920 18.70 1368.00 1.0000000 -1910 26.00 1368.00 1.0000000 -1900 33.80 1368.00 1.0000000 -1890 29.50 1368.00 1.0000000 -1880 18.20 1368.00 1.0000000 -1870 14.40 1368.00 1.0000000 -1860 20.00 1368.00 1.0000000 -1850 26.90 1368.00 1.0000000 -1840 29.00 1368.00 1.0000000 -1830 29.50 1368.00 1.0000000 -1820 34.10 1368.00 1.0000000 -1810 39.70 1368.00 1.0000000 -1800 38.30 1368.00 1.0000000 -1790 33.50 1368.00 1.0000000 -1780 29.40 1368.00 1.0000000 -1770 23.50 1368.00 1.0000000 -1760 17.40 1368.00 1.0000000 -1750 12.60 1368.00 1.0000000 -1740 14.80 1368.00 1.0000000 -1730 25.90 1368.00 1.0000000 -1720 35.30 1368.00 1.0000000 -1710 42.70 1368.00 1.0000000 -1700 47.10 1368.00 1.0000000 -1690 37.70 1368.00 1.0000000 -1680 28.60 1368.00 1.0000000 -1670 25.40 1368.00 1.0000000 -1660 24.50 1368.00 1.0000000 -1650 28.50 1368.00 1.0000000 -1640 34.60 1368.00 1.0000000 -1630 34.40 1368.00 1.0000000 -1620 29.40 1368.00 1.0000000 -1610 27.00 1368.00 1.0000000 -1600 24.60 1368.00 1.0000000 -1590 18.90 1368.00 1.0000000 -1580 12.90 1368.00 1.0000000 -1570 11.50 1368.00 1.0000000 -1560 16.60 1368.00 1.0000000 -1550 21.70 1368.00 1.0000000 -1540 23.20 1368.00 1.0000000 -1530 26.40 1368.00 1.0000000 -1520 31.40 1368.00 1.0000000 -1510 34.40 1368.00 1.0000000 -1500 37.40 1368.00 1.0000000 -1490 39.50 1368.00 1.0000000 -1480 32.50 1368.00 1.0000000 -1470 21.40 1368.00 1.0000000 -1460 18.90 1368.00 1.0000000 -1450 21.50 1368.00 1.0000000 -1440 23.60 1368.00 1.0000000 -1430 27.50 1368.00 1.0000000 -1420 29.70 1368.00 1.0000000 -1410 25.50 1368.00 1.0000000 -1400 21.50 1368.00 1.0000000 -1390 22.00 1368.00 1.0000000 -1380 23.50 1368.00 1.0000000 -1370 21.10 1368.00 1.0000000 -1360 14.40 1368.00 1.0000000 -1350 8.00 1368.00 1.0000000 -1340 5.60 1368.00 1.0000000 -1330 5.60 1368.00 1.0000000 -1320 5.10 1368.00 1.0000000 -1310 7.40 1368.00 1.0000000 -1300 11.10 1368.00 1.0000000 -1290 10.40 1368.00 1.0000000 -1280 11.00 1368.00 1.0000000 -1270 18.70 1368.00 1.0000000 -1260 28.70 1368.00 1.0000000 -1250 29.20 1368.00 1.0000000 -1240 18.10 1368.00 1.0000000 -1230 8.90 1368.00 1.0000000 -1220 10.50 1368.00 1.0000000 -1210 17.10 1368.00 1.0000000 -1200 19.10 1368.00 1.0000000 -1190 19.40 1368.00 1.0000000 -1180 21.70 1368.00 1.0000000 -1170 26.20 1368.00 1.0000000 -1160 30.40 1368.00 1.0000000 -1150 30.50 1368.00 1.0000000 -1140 29.30 1368.00 1.0000000 -1130 25.80 1368.00 1.0000000 -1120 16.90 1368.00 1.0000000 -1110 9.70 1368.00 1.0000000 -1100 12.80 1368.00 1.0000000 -1090 25.10 1368.00 1.0000000 -1080 34.10 1368.00 1.0000000 -1070 32.40 1368.00 1.0000000 -1060 29.20 1368.00 1.0000000 -1050 27.30 1368.00 1.0000000 -1040 24.70 1368.00 1.0000000 -1030 20.60 1368.00 1.0000000 -1020 16.20 1368.00 1.0000000 -1010 14.90 1368.00 1.0000000 -1000 14.70 1368.00 1.0000000 -990 10.30 1368.00 1.0000000 -980 4.40 1368.00 1.0000000 -970 2.70 1368.00 1.0000000 -960 4.70 1368.00 1.0000000 -950 7.90 1368.00 1.0000000 -940 10.60 1368.00 1.0000000 -930 15.70 1368.00 1.0000000 -920 24.40 1368.00 1.0000000 -910 31.00 1368.00 1.0000000 -900 30.30 1368.00 1.0000000 -890 24.60 1368.00 1.0000000 -880 25.10 1368.00 1.0000000 -870 35.80 1368.00 1.0000000 -860 36.10 1368.00 1.0000000 -850 21.70 1368.00 1.0000000 -840 17.70 1368.00 1.0000000 -830 21.50 1368.00 1.0000000 -820 22.40 1368.00 1.0000000 -810 23.10 1368.00 1.0000000 -800 21.40 1368.00 1.0000000 -790 16.20 1368.00 1.0000000 -780 16.10 1368.00 1.0000000 -770 19.90 1368.00 1.0000000 -760 20.90 1368.00 1.0000000 -750 18.80 1368.00 1.0000000 -740 15.20 1368.00 1.0000000 -730 9.30 1368.00 1.0000000 -720 2.10 1368.00 1.0000000 -710 -3.00 1368.00 1.0000000 -700 -3.40 1368.00 1.0000000 -690 -1.50 1368.00 1.0000000 -680 -0.40 1368.00 1.0000000 -670 5.10 1368.00 1.0000000 -660 16.00 1368.00 1.0000000 -650 25.30 1368.00 1.0000000 -640 30.90 1368.00 1.0000000 -630 27.70 1368.00 1.0000000 -620 15.80 1368.00 1.0000000 -610 7.00 1368.00 1.0000000 -600 3.20 1368.00 1.0000000 -590 1.90 1368.00 1.0000000 -580 0.00 1368.00 1.0000000 -570 -3.80 1368.00 1.0000000 -560 -5.50 1368.00 1.0000000 -550 -3.20 1368.00 1.0000000 -540 -0.30 1368.00 1.0000000 -530 -0.30 1368.00 1.0000000 -520 -0.80 1368.00 1.0000000 -510 1.50 1368.00 1.0000000 -500 3.70 1368.00 1.0000000 -490 2.30 1368.00 1.0000000 -480 0.80 1368.00 1.0000000 -470 3.70 1368.00 1.0000000 -460 9.50 1368.00 1.0000000 -450 14.40 1368.00 1.0000000 -440 17.60 1368.00 1.0000000 -430 20.20 1368.00 1.0000000 -420 24.80 1368.00 1.0000000 -410 30.60 1368.00 1.0000000 -400 30.80 1368.00 1.0000000 -390 24.10 1368.00 1.0000000 -380 16.60 1368.00 1.0000000 -370 9.80 1368.00 1.0000000 -360 4.80 1368.00 1.0000000 -350 1.40 1368.00 1.0000000 -340 -2.30 1368.00 1.0000000 -330 -5.10 1368.00 1.0000000 -320 -5.60 1368.00 1.0000000 -310 -4.20 1368.00 1.0000000 -300 -0.70 1368.00 1.0000000 -290 6.70 1368.00 1.0000000 -280 17.80 1368.00 1.0000000 -270 26.10 1368.00 1.0000000 -260 27.30 1368.00 1.0000000 -250 28.30 1368.00 1.0000000 -240 34.40 1368.00 1.0000000 -230 45.30 1368.00 1.0000000 -220 45.90 1368.00 1.0000000 -210 28.90 1368.00 1.0000000 -200 15.20 1368.00 1.0000000 -190 14.00 1368.00 1.0000000 -180 21.50 1368.00 1.0000000 -170 30.50 1368.00 1.0000000 -160 37.50 1368.00 1.0000000 -150 48.70 1368.00 1.0000000 -140 55.00 1368.00 1.0000000 -130 40.00 1368.00 1.0000000 -120 27.00 1368.00 1.0000000 -110 28.70 1368.00 1.0000000 -100 0.00 1368.00 1.0000000 -90 0.00 1368.00 1.0000000 -80 0.00 1368.00 1.0000000 -70 0.00 1368.00 1.0000000 -60 0.00 1368.00 1.0000000 -50 0.00 1368.00 1.0000000 -40 0.00 1368.00 1.0000000 -30 0.00 1368.00 1.0000000 -20 0.00 1368.00 1.0000000 -10 0.00 1368.00 1.0000000 0 0.00 1368.00 1.0000000 10 0.00 1368.00 1.0000000 20 0.00 1368.00 1.0000000 30 0.00 1368.00 1.0000000 40 0.00 1368.00 1.0000000 50 0.00 1368.00 1.0000000 60 0.00 1368.00 1.0000000 70 0.00 1368.00 1.0000000 80 0.00 1368.00 1.0000000 90 0.00 1368.00 1.0000000 100 0.00 1368.00 1.0000000 110 0.00 1368.00 1.0000000 120 0.00 1368.00 1.0000000 130 0.00 1368.00 1.0000000 140 0.00 1368.00 1.0000000 150 0.00 1368.00 1.0000000 160 0.00 1368.00 1.0000000 170 0.00 1368.00 1.0000000 180 0.00 1368.00 1.0000000 190 0.00 1368.00 1.0000000 200 0.00 1368.00 1.0000000 210 0.00 1368.00 1.0000000 220 0.00 1368.00 1.0000000 230 0.00 1368.00 1.0000000 240 0.00 1368.00 1.0000000 250 0.00 1368.00 1.0000000 260 0.00 1368.00 1.0000000 270 0.00 1368.00 1.0000000 280 0.00 1368.00 1.0000000 290 0.00 1368.00 1.0000000 300 0.00 1368.00 1.0000000 310 0.00 1368.00 1.0000000 320 0.00 1368.00 1.0000000 330 0.00 1368.00 1.0000000 340 0.00 1368.00 1.0000000 350 0.00 1368.00 1.0000000 360 0.00 1368.00 1.0000000 370 0.00 1368.00 1.0000000 380 0.00 1368.00 1.0000000 390 0.00 1368.00 1.0000000 400 0.00 1368.00 1.0000000 410 0.00 1368.00 1.0000000 420 0.00 1368.00 1.0000000 430 0.00 1368.00 1.0000000 440 0.00 1368.00 1.0000000 450 0.00 1368.00 1.0000000 460 0.00 1368.00 1.0000000 470 0.00 1368.00 1.0000000 480 0.00 1368.00 1.0000000 490 0.00 1368.00 1.0000000 500 0.00 1368.00 1.0000000 510 0.00 1368.00 1.0000000 520 0.00 1368.00 1.0000000 530 0.00 1368.00 1.0000000 540 0.00 1368.00 1.0000000 550 0.00 1368.00 1.0000000 560 0.00 1368.00 1.0000000 570 0.00 1368.00 1.0000000 580 0.00 1368.00 1.0000000 590 0.00 1368.00 1.0000000 600 0.00 1368.00 1.0000000 610 0.00 1368.00 1.0000000 620 0.00 1368.00 1.0000000 630 0.00 1368.00 1.0000000 640 0.00 1368.00 1.0000000 650 0.00 1368.00 1.0000000 660 0.00 1368.00 1.0000000 670 0.00 1368.00 1.0000000 680 0.00 1368.00 1.0000000 690 0.00 1368.00 1.0000000 700 0.00 1368.00 1.0000000 710 0.00 1368.00 1.0000000 720 0.00 1368.00 1.0000000 730 0.00 1368.00 1.0000000 740 0.00 1368.00 1.0000000 750 0.00 1368.00 1.0000000 760 0.00 1368.00 1.0000000 770 0.00 1368.00 1.0000000 780 0.00 1368.00 1.0000000 790 0.00 1368.00 1.0000000 800 0.00 1368.00 1.0000000 810 0.00 1368.00 1.0000000 820 0.00 1368.00 1.0000000 830 0.00 1368.00 1.0000000 840 0.00 1368.00 1.0000000 850 0.00 1368.00 1.0000000 860 0.00 1368.00 1.0000000 870 0.00 1368.00 1.0000000 880 0.00 1368.00 1.0000000 890 0.00 1368.00 1.0000000 900 0.00 1368.00 1.0000000 910 0.00 1368.00 1.0000000 920 0.00 1368.00 1.0000000 930 0.00 1368.00 1.0000000 940 0.00 1368.00 1.0000000 950 0.00 1368.00 1.0000000 960 0.00 1368.00 1.0000000 970 0.00 1368.00 1.0000000 980 0.00 1368.00 1.0000000 990 0.00 1368.00 1.0000000 1000 0.00 1368.00 1.0000000 1010 0.00 1368.00 1.0000000 1020 0.00 1368.00 1.0000000 1030 0.00 1368.00 1.0000000 1040 0.00 1368.00 1.0000000 1050 0.00 1368.00 1.0000000 1060 0.00 1368.00 1.0000000 1070 0.00 1368.00 1.0000000 1080 0.00 1368.00 1.0000000 1090 0.00 1368.00 1.0000000 1100 0.00 1368.00 1.0000000 1110 0.00 1368.00 1.0000000 1120 0.00 1368.00 1.0000000 1130 0.00 1368.00 1.0000000 1140 0.00 1368.00 1.0000000 1150 0.00 1368.00 1.0000000 1160 0.00 1368.00 1.0000000 1170 0.00 1368.00 1.0000000 1180 0.00 1368.00 1.0000000 1190 0.00 1368.00 1.0000000 1200 0.00 1368.00 1.0000000 1210 0.00 1368.00 1.0000000 1220 0.00 1368.00 1.0000000 1230 0.00 1368.00 1.0000000 1240 0.00 1368.00 1.0000000 1250 0.00 1368.00 1.0000000 1260 0.00 1368.00 1.0000000 1270 0.00 1368.00 1.0000000 1280 0.00 1368.00 1.0000000 1290 0.00 1368.00 1.0000000 1300 0.00 1368.00 1.0000000 1310 0.00 1368.00 1.0000000 1320 0.00 1368.00 1.0000000 1330 0.00 1368.00 1.0000000 1340 0.00 1368.00 1.0000000 1350 0.00 1368.00 1.0000000 1360 0.00 1368.00 1.0000000 1370 0.00 1368.00 1.0000000 1380 0.00 1368.00 1.0000000 1390 0.00 1368.00 1.0000000 1400 0.00 1368.00 1.0000000 1410 0.00 1368.00 1.0000000 1420 0.00 1368.00 1.0000000 1430 0.00 1368.00 1.0000000 1440 0.00 1368.00 1.0000000 1450 0.00 1368.00 1.0000000 1460 0.00 1368.00 1.0000000 1470 0.00 1368.00 1.0000000 1480 0.00 1368.00 1.0000000 1490 0.00 1368.00 1.0000000 1500 0.00 1368.00 1.0000000 1510 0.00 1368.00 1.0000000 1520 0.00 1368.00 1.0000000 1530 0.00 1368.00 1.0000000 1540 0.00 1368.00 1.0000000 1550 0.00 1368.00 1.0000000 1560 0.00 1368.00 1.0000000 1570 0.00 1368.00 1.0000000 1580 0.00 1368.00 1.0000000 1590 0.00 1368.00 1.0000000 1600 0.00 1368.00 1.0000000 1610 0.00 1368.00 1.0000000 1620 0.00 1368.00 1.0000000 1630 0.00 1368.00 1.0000000 1640 0.00 1368.00 1.0000000 1650 0.00 1368.00 1.0000000 1660 0.00 1368.00 1.0000000 1670 0.00 1368.00 1.0000000 1680 0.00 1368.00 1.0000000 1690 0.00 1368.00 1.0000000 1700 0.00 1368.00 1.0000000 1710 0.00 1368.00 1.0000000 1720 0.00 1368.00 1.0000000 1730 0.00 1368.00 1.0000000 1740 0.00 1368.00 1.0000000 1750 0.00 1368.00 1.0000000 1760 0.00 1368.00 1.0000000 1770 0.00 1368.00 1.0000000 1780 0.00 1368.00 1.0000000 1790 0.00 1368.00 1.0000000 1800 0.00 1368.00 1.0000000 1810 0.00 1368.00 1.0000000 1820 0.00 1368.00 1.0000000 1830 0.00 1368.00 1.0000000 1840 0.00 1368.00 1.0000000 1850 0.00 1368.00 1.0000000 1860 0.00 1368.00 1.0000000 1870 0.00 1368.00 1.0000000 1880 0.00 1368.00 1.0000000 1890 0.00 1368.00 1.0000000 1900 0.00 1368.00 1.0000000 1910 0.00 1368.00 1.0000000 1920 0.00 1368.00 1.0000000 1930 0.00 1368.00 1.0000000 1940 0.00 1368.00 1.0000000 1950 0.00 1368.00 1.0000000 1960 0.00 1368.00 1.0000000 1970 0.00 1368.00 1.0000000 1980 0.00 1368.00 1.0000000 1990 0.00 1368.00 1.0000000 2000 0.00 1368.00 1.0000000 2010 0.00 1368.00 1.0000000 2020 0.00 1368.00 1.0000000 2030 0.00 1368.00 1.0000000 2040 0.00 1368.00 1.0000000 2050 0.00 1368.00 1.0000000 2060 0.00 1368.00 1.0000000 2070 0.00 1368.00 1.0000000 2080 0.00 1368.00 1.0000000 2090 0.00 1368.00 1.0000000 2100 0.00 1368.00 1.0000000 2110 0.00 1368.00 1.0000000 2120 0.00 1368.00 1.0000000 2130 0.00 1368.00 1.0000000 2140 0.00 1368.00 1.0000000 2150 0.00 1368.00 1.0000000 2160 0.00 1368.00 1.0000000 2170 0.00 1368.00 1.0000000 2180 0.00 1368.00 1.0000000 2190 0.00 1368.00 1.0000000 2200 0.00 1368.00 1.0000000 2210 0.00 1368.00 1.0000000 2220 0.00 1368.00 1.0000000 2230 0.00 1368.00 1.0000000 2240 0.00 1368.00 1.0000000 2250 0.00 1368.00 1.0000000 2260 0.00 1368.00 1.0000000 2270 0.00 1368.00 1.0000000 2280 0.00 1368.00 1.0000000 2290 0.00 1368.00 1.0000000 2300 0.00 1368.00 1.0000000 2310 0.00 1368.00 1.0000000 2320 0.00 1368.00 1.0000000 2330 0.00 1368.00 1.0000000 2340 0.00 1368.00 1.0000000 2350 0.00 1368.00 1.0000000 2360 0.00 1368.00 1.0000000 2370 0.00 1368.00 1.0000000 2380 0.00 1368.00 1.0000000 2390 0.00 1368.00 1.0000000 2400 0.00 1368.00 1.0000000 2410 0.00 1368.00 1.0000000 2420 0.00 1368.00 1.0000000 2430 0.00 1368.00 1.0000000 2440 0.00 1368.00 1.0000000 2450 0.00 1368.00 1.0000000 2460 0.00 1368.00 1.0000000 2470 0.00 1368.00 1.0000000 2480 0.00 1368.00 1.0000000 2490 0.00 1368.00 1.0000000 2500 0.00 1368.00 1.0000000 2510 0.00 1368.00 1.0000000 2520 0.00 1368.00 1.0000000 2530 0.00 1368.00 1.0000000 2540 0.00 1368.00 1.0000000 2550 0.00 1368.00 1.0000000 2560 0.00 1368.00 1.0000000 2570 0.00 1368.00 1.0000000 2580 0.00 1368.00 1.0000000 2590 0.00 1368.00 1.0000000 2600 0.00 1368.00 1.0000000 ; expanding heap to 240000 words Solanki_rows := first floor ((shape Solanki) / Solanki_cols ) ; Solanki := Solanki_rows Solanki_cols reshape Solanki ; # should double-check to make sure "Solanki_rows Solanki_cols" have been changed if data changes # enddoc # Civilization problem setup.ndf 10Feb07 www.BillHowell.ca # ********************************* # Original code from: J. Laskar, P. Robutel, F. Joutel, M. Gastineau, A.C.M. Correia, B. Levrard 2004 "A long term numerical solution for the insolation quantities of the Earth" Astronomy & Astrophysics manuscript no. La 2004 May 23, 2004 http://www.imcce.fr/fr/presentation/equipes/ASD/preprints/prep.2004/La_2004_prep.pdf # Laskar, J., Joutel, F., Boudin, F. 1993 "Orbital, precessional and insolation quantities for the Earth from -20 Myr to + 10Myr" Astron. Astrophys. 270, 522 http://www.imcce.fr/Equipes/ASD/insola/earth/La93/README.TXT # Earth coordinates read APPROXIMATELY from MapQuest: http://www.mapquest.com/atlas/ Quick background on civilizations via www.Wikipedia.com (great source of information) # Warning: Use Courier constant-width font and no-line-wrap to view this file! Otherwise the code doesn't line up and it's hard to read. # Just for a standard range of dates for all civilizations (past 11 ky out to 2 ky into future, but include ~10 points either side to ensure polynomial 8th order interpolation is good). Note that Laskar etal define start & end times in terms of millions of years but it is more "human" to use years for the extremely short timeframes here plus this avoids some confusion. These variables are set (checked) in the Milankovic parameter file: datedebut % starting time (years) ; datefin % final time (years) ; pas % stepsize (years) must be greater than 1 year ; # All of the insolation calculations below are specified for 21Jun-21Jul insolation ; More examples with greater flexibility will be added later in different parameter files WARNING: remove apostrophes from the longitude/latitude coordinates!!! # approximate conventional dates for a given mean longitude' 0: 21 march ' 30: 21 april ' 60: 21 may ' 90: 21 june ' 120: 21 july ' 150: 21 august ' 180: 21 september' 210: 21 october ' 240: 21 november ' 270: 21 december ' 300: 21 january ' 330: 21 february ' # approximate conventional dates of the months' 1: 21 december - 20 january' 2: 21 january - 20 february' 3: 21 february - 20 march' 4: 21 march - 20 april' 5: 21 april - 20 may' 6: 21 may - 20 june' 7: 21 june - 20 july' 8: 21 july - 20 august' 9: 21 august - 20 september' 10: 21 september - 20 october' 11: 21 october - 20 november' 12: 21 november - 20 december' # ********************************* # globals defined here, or used in this file only ; civic_e := f99 ; # EXTERNAL definitions to allow this file to load alone civilization_insolation IS external operation ; civilization_insolation_x := f99 ; complete_inputs IS external expression ; complete_inputs_x := f99 ; # ************************** ; % Check that all user inputs are OK, and calculate intermediate variables ; Civic IS { IF complete_inputs THEN % Mayan civilization calcs ; fichout := link Laskar_lib 'Mayan insolation Jun-Jul.txt' ; % MapQuest ~21N 89.5W Mérida on Yucatan peninsula ; latitude := 21 / 180 ; % Mayan latitude in radians North - positive, South - negative latitude ; mois := 7 ; % mean monthly, 7 => 21 june - 20 july ; civilization_insolation fichout latitude mois ; % Wikipedia nyet ; % Roman empire calcs ; % Wikipedia ; % MapQuest ~42N 12.5E Rome ; % ; % Anasazi indians in SW USA ; % Wikipedia Ancient Pueblo Peoples http://en.wikipedia.org/wiki/Anasazi ; % four U.S. states of Utah, Colorado, New Mexico, and Arizona ; % MapQuest ~35N 109W in middle of 4 states ; % ; % Harrapa civilization in the Indus valley (!&!Pakistan!&!) ; % Wikipedia ; % MapQuest ~25N 68.5E Hyderabad ; % ; % Old Kingdom in Egypt ; % Wikipedia - The Old Kingdom is most commonly regarded as spanning the period ; % of time when Egypt was ruled by the Third Dynasty through to the ; % Sixth Dynasty (2686 BC - 2134 BC). Capital city Memphis ; % MapQuest ~30N 31E Memphis ; % ; % Easter Island ; % Wikipedia It is located at 27°09S 109°27Wt is one of the most ; % isolated inhabited islands in the world ; % ; % Minoan civilization ; % Wikipedia puts on island of Crete, 2700 to 1450 BC ; % MapQuest ~35N 25E ; % ; % Mesopotamia - Sumeria, Akkadia, Babylonian ; % MapQuest ~31N 47.5E Tigis & Euphrates juncture ; % ; % China - !&!Tien state!&! ; % MapQuest ~40N 116E Beijing ; % ; % ; ELSE writescreen 'Civilization problem setup.ndf: Error returned by complete_inputs!!' ; ENDIF ; } # enddoc # Main.ndf 17Feb07 www.BillHowell.ca Warning: Use Courier constant-width font and no-line-wrap to view this file! Otherwise the code doesn't line up and it's hard to read. # Original code from: J. Laskar, P. Robutel, F. Joutel, M. Gastineau, A.C.M. Correia, B. Levrard 2004 "A long term numerical solution for the insolation quantities of the Earth" Astronomy & Astrophysics manuscript no. La 2004 May 23, 2004 http://www.imcce.fr/fr/presentation/equipes/ASD/preprints/prep.2004/La_2004_prep.pdf # Laskar, J., Joutel, F., Boudin, F. 1993 "Orbital, precessional and insolation quantities for the Earth from -20 Myr to + 10Myr" Astron. Astrophys. 270, 522 http://www.imcce.fr/Equipes/ASD/insola/earth/La93/README.TXT #**************************************** * * traduction en Q'Nial de www.qnial.com - * Michael, Jenkins & team at Queen's University, Canada * par www.BillHowell.ca février 2007 * * WARNING: This a draft port from fortran to Q'Nial. Errors are guaranteed * * USER INPUTS: Rather than use a binary file for database access to intermediate results, * I've simply used an array in memory (my current focus is limited to the Holocene period). * I'm not sure if an array of 100 k records of 5 to 10 reals each would be too * much memory for this version of Q'Nial, but I think it would be OK. * * Load definitions - with Q'Nial, it is the problem spec file that will drive the actions * (Howell - I prefer not to waste my time with user input prompts, although that option * will eventually be supported) * * Insolation operators (cwj, wmcal, wjour, wjcal, Global_annual_avg_insolation) return a real, * they don't assign to an argument as in theoriginal Laskar etal * * Missing (some other time): * help files with search - in pdf format and pretty fonts etc etc (yuuchh, but that's what people want) * wide range of examples for different types of research * * * Other comments for Q'Nial speedup * replace "x power 2" with "x * x" * don't keep recreating arrays !!! -> QROMB POLINT * * * Bill Howell * #********************************* # globals defined here, or used in this file only ; datedebut := 0.0 ; datedebut_v := f1 ; datefin := 0.0 ; datefin_v := f1 ; itdebut := 0.0 ; itdebut_v := f1 ; % Nombre de lignes lues ; aux := 0.0 ; aux_v := f1 ; npt := 0.0 ; npt_v := f1 ; % Milankovic_rows !! ; civilisation_period := 0.0 ; civilisation_period_v := f1 ; Milankovic_delta_t := 0.0 ; Milankovic_delta_t_v := f1 ; Solanki_delta_t := 0.0 ; Solanki_delta_t_v := f1 ; Complete_inputs_e := f1 ; # EXTERNAL definitions to allow this file to load alone date_conversion IS external variable ; date_conversion_x := f1 ; debut IS external variable ; debut_x := f1 ; interpolation_setup IS external expression ; interpolation_setup_x := f1 ; Milankovic IS external variable ; Milankovic_x := f1 ; Milankovic_rows IS external variable ; Milankovic_rows_x := f1 ; n_intrpsteps IS external variable ; n_intrpsteps_x := f1 ; #**************************************** # Set up files, vérification du pas Complete_inputs IS { NONLOCAL civilisation_period datedebut datefin itdebut npt ; % ; % Milankovic columns ; yr_col := 0 ; ecc_col := 1 ; obl_col := 2 ; cosp_col := 3 ; sinp_col := 4 ; prec_col := 5 ; input_error := l ; % Extra care to make user-define variables are all there!!! ; symbol_list := EACH first symbols 0 ; IF not in "MILANKOVIC symbol_list THEN writescreen '?Error in file: Civilization Milankovic precess, obliq, eccent.txt: ' ; writescreen 'Milankovic data hasnt been loaded' ; input_error := l ; ENDIF ; IF not in "DATE_CONVERSION symbol_list THEN writescreen '?Error in file: Civilization Milankovic precess, obliq, eccent.txt: ' ; writescreen 'Missing date_coversion to change first column of data input to years from ky or My' ; input_error := l ; ENDIF ; IF not in "PAS symbol_list THEN writescreen '?Error in file: Civilization Milankovic precess, obliq, eccent.txt: ' ; writescreen 'Missing stepsize (years - NOT millions of years)' ; input_error := l ; ENDIF ; % data completion & testing continued ; Milankovic|[, yr_col ] := Milankovic|[ , yr_col ] * date_conversion + 2000 ; Milankovic|[, cosp_col] := cos Milankovic|[ , prec_col] ; Milankovic|[, sinp_col] := sin Milankovic|[ , prec_col] ; datedebut := Milankovic@[0, yr_col ] ; % starting time (years) ; datefin := Milankovic@[Milankovic_rows - 1,0] ; % final time (years) ; % Fancy data checks ; % There must be constant time intervals between data dates!! ; IF ~= (Milankovic_delta_t := - [rest, front] Milankovic|[,0]) THEN writescreen '?Error in file: Milankovic - civilization.ndf' ; writescreen 'Milankovic data is provided for an irregular date series - there must be constant years difference between data in the series' ; input_error := l ; ELSE Milankovic_delta_t := first Milankovic_delta_t ; ENDIF ; IF Milankovic_delta_t < 1 THEN writescreen '?Error in file: Milankovic - civilization.ndf ' ; writescreen 'Step size (Milankovic_delta_t) must be >= 1 year' ; input_error := l ; ENDIF ; IF NOT < (datedebut datefin) THEN writescreen '?Error in file: Civilization Milankovic precess, obliq, eccent.txt: ' ; writescreen 'Milankovic dates must start with earliest and end with last' ; input_error := l ; ENDIF ; IF ~= (Solanki_delta_t := - [rest, front] Solanki|[,0]) THEN writescreen '?Error in file: Solanki-Tapping.ndf' ; writescreen 'Solanki-Tapping data is provided for an irregular date series - there must be constant years difference between data in the series' ; input_error := l ; ELSE Solanki_delta_t := first Solanki_delta_t ; ENDIF ; interpolation_setup ; % I can't think of any checks after initial program debugging ; % Complete the input data by doing extra calculations ; civilisation_period := -11000 + (1000 * tell 12) ; NOT input_error }