attachment:darkside0.pl of Brown2013

   1 #!/usr/local/bin/perl
   2 # darkside0.pl
   3 
   4 use Math::Trig;
   5 
   6 my $day     = 23*3600 + 56*60 + 4.090530833  ; # seconds per sidereal day
   7 
   8 my $comppts = 100           ; # computations per plot
   9 my $plotpts = 500           ; # plot points
  10 my $tmax    = 5000          ; # max seconds
  11 
  12 my $m0      = 24*60         ; # minutes per synodic day
  13 my $mu      = 398600.4418   ; # km^3/s^2
  14 my $pi2     = 8.0*atan(1.0) ; # 2 pi
  15 my $re      = 6378.1        ; # km
  16 my $etemp   = 250.0         ; # Earth temp in Kelvin
  17 my $stblz   = 5.6704e-08    ; # W/m2-K4 Stephan Boltzmann constant
  18 my $sun     = 1300.0        ; # Sun watts/m2
  19 my $mass    = 0.300         ; # kg/m2
  20 my $spheat  = 710.0         ; # Specific heat J/kg-K
  21 my $albedo  = 0.85          ; # serversat albedo
  22 my $celsius = -273.15       ; # kelvin to celsius
  23 
  24 my $deltat  = $tmax / ($plotpts*$comppts ) ;
  25 my $heatcap = $mass * $spheat  ;
  26 my $eshine  = $stblz*$etemp**4 ; # earth shine in watts/
  27 
  28 my @array ;
  29 
  30 for( my $n=1 ; $n < 6 ; $n++ ) {
  31    my $norb       = 6 - $n       ;
  32    my $nper       = 7 - $n       ;
  33    $nper = 1      if $nper == 2  ; 
  34    my $min        = $m0  / $norb ;
  35    my $period     = $day / $nper ;
  36    my $semimajor  = ( $mu * ( $period / $pi2 )**2.0 )**(1.0/3.0) ;
  37    my $darkangle  = 2.0 * asin( $re / $semimajor ) ;
  38    my $darktime   = $period * $darkangle / $pi2 ;
  39    my $earthsky   = 0.5 * ( 1.0 - cos( 0.5*$darkangle ) ) ;
  40    
  41    # 2 surfaces
  42    my $earthwatts = $albedo * $earthsky*$stblz*$etemp**4 ; 
  43    my $startwatts = $albedo * $sun + $earthwatts ;
  44    my $starttemp  = sqrt(sqrt($startwatts/(2.0*$albedo*$stblz) ) ) ;
  45 
  46    ## printf "m%04d%9.3f%9.1f%9.5f%9.2f\n",
  47    ##          $min, $starttemp, $semimajor, $darkangle, $darktime ;
  48 
  49    my $time  = 0           ;
  50    my $temp  = $starttemp  ;
  51    my $watts = $earthwatts ;
  52       $array[ 0][0] = $time ;
  53       $array[$n][0] = $temp + $celsius ;
  54    for( my $pt = 1 ; $pt < $plotpts ; $pt++ ) {
  55       for( my $comp = 0 ; $comp < $comppts ; $comp++ ) {
  56          my $wattsin = $watts - 2.0*$albedo*$stblz*$temp**4 ;
  57          $temp += $wattsin * $deltat / $heatcap ;
  58          $time += $deltat ;
  59          $watts = $startwatts if ( $time > $darktime ) ;
  60       }
  61       $array[ 0][$pt] = $time ;
  62       $array[$n][$pt] = $temp + $celsius ;
  63 }  }  
  64 
  65 
  66 for( my $pt = 0 ; $pt < $plotpts ; $pt++ ) {
  67    for( my $n=0 ; $n < 6 ; $n++ ) {
  68       printf "%9.3f", $array[$n][$pt] ;
  69    }
  70    printf "\n" ;
  71 }