#!/usr/local/bin/perl
# darkside0.pl

use Math::Trig;

my $day     = 23*3600 + 56*60 + 4.090530833  ; # seconds per sidereal day

my $comppts = 100           ; # computations per plot
my $plotpts = 500           ; # plot points
my $tmax    = 5000          ; # max seconds

my $m0      = 24*60         ; # minutes per synodic day
my $mu      = 398600.4418   ; # km^3/s^2
my $pi2     = 8.0*atan(1.0) ; # 2 pi
my $re      = 6378.1        ; # km
my $etemp   = 250.0         ; # Earth temp in Kelvin
my $stblz   = 5.6704e-08    ; # W/m2-K4 Stephan Boltzmann constant
my $sun     = 1300.0        ; # Sun watts/m2
my $mass    = 0.300         ; # kg/m2
my $spheat  = 710.0         ; # Specific heat J/kg-K
my $albedo  = 0.85          ; # serversat albedo
my $celsius = -273.15       ; # kelvin to celsius

my $deltat  = $tmax / ($plotpts*$comppts ) ;
my $heatcap = $mass * $spheat  ;
my $eshine  = $stblz*$etemp**4 ; # earth shine in watts/

my @array ;

for( my $n=1 ; $n < 6 ; $n++ ) {
   my $norb       = 6 - $n       ;
   my $nper       = 7 - $n       ;
   $nper = 1      if $nper == 2  ; 
   my $min        = $m0  / $norb ;
   my $period     = $day / $nper ;
   my $semimajor  = ( $mu * ( $period / $pi2 )**2.0 )**(1.0/3.0) ;
   my $darkangle  = 2.0 * asin( $re / $semimajor ) ;
   my $darktime   = $period * $darkangle / $pi2 ;
   my $earthsky   = 0.5 * ( 1.0 - cos( 0.5*$darkangle ) ) ;
   
   # 2 surfaces
   my $earthwatts = $albedo * $earthsky*$stblz*$etemp**4 ; 
   my $startwatts = $albedo * $sun + $earthwatts ;
   my $starttemp  = sqrt(sqrt($startwatts/(2.0*$albedo*$stblz) ) ) ;

   ## printf "m%04d%9.3f%9.1f%9.5f%9.2f\n",
   ##          $min, $starttemp, $semimajor, $darkangle, $darktime ;

   my $time  = 0           ;
   my $temp  = $starttemp  ;
   my $watts = $earthwatts ;
      $array[ 0][0] = $time ;
      $array[$n][0] = $temp + $celsius ;
   for( my $pt = 1 ; $pt < $plotpts ; $pt++ ) {
      for( my $comp = 0 ; $comp < $comppts ; $comp++ ) {
         my $wattsin = $watts - 2.0*$albedo*$stblz*$temp**4 ;
         $temp += $wattsin * $deltat / $heatcap ;
         $time += $deltat ;
         $watts = $startwatts if ( $time > $darktime ) ;
      }
      $array[ 0][$pt] = $time ;
      $array[$n][$pt] = $temp + $celsius ;
}  }  


for( my $pt = 0 ; $pt < $plotpts ; $pt++ ) {
   for( my $n=0 ; $n < 6 ; $n++ ) {
      printf "%9.3f", $array[$n][$pt] ;
   }
   printf "\n" ;
}