### dw.tcl, cr sun 07 may 2000 by rha
# to simulate alpha vs L

### set our process constants##################
# this one prints debug statements
set DEBUG 0
# this one prints data for animation routine
set ANIM 1

### set constants from WL ###
## these are GLOBAL constants

# flux constant
set S 9.17e+5	
# set stefan's constant
# set sigma 0.000056 
set sigma 0.0000599 

## these are GLOBAL parameters

# death rate both daisies
set gamma 0.3	
# total fertile area
set P 1.0	
# mixing coefficient
set qprime 20.0	

# test them
if {$DEBUG} {
  puts "WL constants (gamma, P, S, sigma, qprime) are $gamma, $P, $S, $sigma, $qprime"
}

# local albedo assumptions
set AB 0.25
set AW 0.75
set AG 0.50
# test them
if {$DEBUG} {
  puts "WL albedos (AB, AW, AG) are $AB, $AW, $AG"
}

# range of L for response diagrams
set Lmin 0.6
set Lmax 1.7
if {$DEBUG} {
  puts "WL range of Lmin, max: $Lmin $Lmax"
}

### inits for first L sim, we guess these initial areas
set alphaBinit 0.6
set alphaWinit 0.2
set xinit [expr {$P - $alphaBinit - $alphaWinit}]
set alphaGinit $xinit
if {$DEBUG} {
  puts "our initial areas are $alphaBinit, $alphaWinit, $alphaGinit"
}

### some funcs from WL #####################

### find albedo, A WL eqn (5)
## assume P=1, let a denote area or alpha
proc albedo {aB aW} {
  global AB AW AG DEBUG
  set aG [expr {1 - $aB - $aW}]
  set alb [expr {($aB * $AB) + ($aW * $AW) + ($aG * $AG)}]
  if {$DEBUG} {
    puts "within proc albedo, alb is, $alb"
  }
  
  return $alb
}

## find effective temperature, Te
# have to solve fourth degree equation, WL(4)
# let TE = Te + 273

proc Te {lum aB aW} {
  global S sigma DEBUG
  set coeff [expr {$S*$lum/$sigma}]
  set A [albedo $aB $aW]
  set myarg [expr {$coeff*(1-$A)}]
  set TE [expr {pow($myarg, 0.25)}]
  set te [expr {$TE - 273.0}]
  # test them
  if {$DEBUG} {
    puts "coeff is $coeff"
    puts "albedo is $A"
    puts "myarg is $myarg"
    puts "eff temp is: $te"
  }
  return $te
}

### the parabolic birth rate function
## from WL eqn (3)

proc beta {t} {
  global DEBUG
  set square [expr { (22.5 - $t)*(22.5 - $t) }]
  set quad [expr {1 - 0.003265*$square}]
  if { $quad > 0 } {
    set b $quad
  } else {
    set b 0.0
  }
  if {$DEBUG} {
    puts "within proc beta, b is: $b"
  }
  return $b  
}

### L loop will begin here ###################
## initial stuff
# fix first L
set m 0
set mmax 11
set deltaL [expr {($Lmax - $Lmin)/$mmax}]
if {$DEBUG} {
  puts "deltaL is: $deltaL"
}

# init areas
set alphaBnow $alphaBinit
set alphaWnow $alphaWinit
set xnow [expr {$P - $alphaBnow - $alphaWnow}]

## now loop
# NOTE: following WL, we do not reinit areas for each m
for {set m 0} {$m < $mmax} {incr m 1} {
	if {$DEBUG} {
	  puts "/n at begin of L loop, areas are: $alphaBnow, $alphaWnow"
	}
	set L [expr {$Lmin + $m * $deltaL}]
	set A [albedo $alphaBnow $alphaWnow]
	set T [Te $L $alphaBnow $alphaWnow]
	## next, local temperatures from linear approx WL(7)
	set TB [expr {$qprime*($A - $AB) + $T}]
	set TW [expr {$qprime*($A - $AW) + $T}]
	if {$DEBUG} {
	  puts "TB, TW are: $TB, $TW"
	}
	set betaB [beta $TB]
	set betaW [beta $TW]

	## finally, the vectorfield!!
	set alphaBrate [expr {$alphaBnow*($xnow * $betaB - $gamma)}]
	set alphaWrate [expr {$alphaWnow*($xnow * $betaW - $gamma)}]
	if {$DEBUG} {
	  puts "alpha rates are now: $alphaBrate, $alphaWrate"
	}

	### euler routine
	## setup euler constants
	# total time, numbersteps, delta time
	set nmax 20
	set dt 0.01
	set tmax [expr {$nmax * $dt}]
	# test them
	if {$DEBUG} {
	  puts "Euler constants are $tmax, $nmax, $dt"
	}

	# now do nmax steps
	for {set i 0} {$i < $nmax} {incr i 1} {
	  set alphaBnext [expr {$alphaBnow + $dt * $alphaBrate}]
	  if {$alphaBnext < 0} {
	    set alphaBnext 0.0
	  }
	  set alphaWnext [expr {$alphaWnow + $dt * $alphaWrate}]
	  if {$alphaWnext < 0} {
	    set alphaWnext 0.0
	  }
	  # puts "after step $i: now, next Bstates are: $alphaBnow, $alphaBnext"
	  # puts "after step $i: now, next Wstates are: $alphaWnow, $alphaWnext"
	  # ready for next step
	  set alphaBnow $alphaBnext
	  set alphaWnow $alphaWnext
	  set xnow [expr {$P - $alphaBnow - $alphaWnow}]
	  set A [albedo $alphaBnow $alphaWnow]
	  set T [Te $L $alphaBnow $alphaWnow]
	  ## next, local temperatures from linear approx WL(7)
	  set TB [expr {$qprime*($A - $AB) + $T}]
	  set TW [expr {$qprime*($A - $AW) + $T}]
	  if {$DEBUG} {
	    puts "TB, TW are: $TB, $TW"
	  }
	  set betaB [beta $TB]
	  set betaW [beta $TW]
	}

	# if ANIM, write to stdout
	if {$ANIM} {
	  puts "$L, $alphaBnow, $alphaWnow, $T"
	}
}

exit

#end: dw.tcl