--- Daisyworld model 

 --[[
 Author:      Gilberto Camara
 Model name:  Daisyworld
 Date:        22.11.2015
 
]]--

--[[

Wood, A. J., G. J. Ackland, J. G. Dyke, H. T. P. Williams, and T. M. Lenton (2008), 
Daisyworld: A review, Rev. Geophys., 46.

Daisyworld consists of two different types of ‘‘daisy,’’ which may be considered distinct species (because there is no possibility of mixed replication of the types) or, alternatively, as distinct phenotypes of the same species. The two types are identified as either black or white according to their reflectivity or albedo.

--]]

--[[
Daisy growth rate function is a function of temperature. Maximum growth 
occurs at 22.5 degrees Celsius, and daisies grow from a minimum of 5 C to 
a maximum of 40 C.
--]]

---------------------------------------------------------------------------------------
function daisyGrowthRate (tempK) -- daisyGrowthRate as function of temperature
    local gr = 0.0
    local temp = tempK - 273
    if temp > 5.0 and temp < 40.0 then 
        gr = 1 - 0.003265 * (22.5 - temp)^2
    end      
    return gr
end

--- function that calculates luminosity 

--[[ 
Calculates the average planet temperature using Boltzmann's law 

"The amount of thermal radiation emitted increases rapidly and the principal frequency of the radiation becomes higher with increasing temperatures" -- 

TempK^4 = average solar flux * luminosity (1 - albedo)

 Watson and Lovelock used 3668 watts per square meter as the solar flux,
 which is much higher than the average solar flux incident on the Earth today, 
 because the Earth has a relatively low planetary albedo and 
 a significant greenhouse effect. 
]]

SIGMA = 5.67E-08    -- Joules / per sec * square meters * TempK^^4 -- Stephen Boltzmann constant
SOLAR_FLUX = 3668   -- Watts/square meter

function calcTemp (lum, albedo)

    local tempK = (lum * SOLAR_FLUX * ( 1 - albedo)/(4*SIGMA)) ^0.25
    return tempK 
end

--[[
Calculates the temperatures near the daisies

The local temperatures Tw and Tb are defined by making a simplifying assumption about 
the heat transfer: a linearization of a diffusion term. 
A parameter q is defined as the heat transfer coefficient, 
thereby defining the local temperatures as

TempK^4 (daisy) = q (planet_albedo - daisy_albedo) + tempK^4 (planet)
--]]

HEAT_TRANSFER_COEF = 2.06E09  -- Kelvins^4  -- used by Watson and Lovelock

function tempNearDaisy (planet_temp, planet_albedo, daisy_albedo)
    local   temp_daisy  = (HEAT_TRANSFER_COEF * (planet_albedo - daisy_albedo) + planet_temp^4 ) ^0.25
    return  temp_daisy
end

--[[
-- actual growthRate is constained by the empty area of the planet 

--]]
function growthReduction (baseRate, emptyArea)
    local newRate =  baseRate * emptyArea
    return newRate
end

--[[ 
    planet albedo
--]]

function calcPlanetAlbedo (model)
    local albedo = model.whiteArea * model.whiteAlbedo + 
                   model.blackArea * model.blackAlbedo +
                   model.emptyArea * model.soilAlbedo 
    return albedo
end


function growDaisies (model)
    
    -- runs the daisyworld model

    -- gets the average albedo of the planet
    model.planetAlbedo = calcPlanetAlbedo (model)
        
    -- function that calculates the planet Temperature
    model.planetTemp = calcTemp (model.sunLuminosity, model.planetAlbedo)
        
    -- temperature near the white daisies
    model.tempNearWhite = tempNearDaisy (model.planetTemp, 
            model.planetAlbedo, model.whiteAlbedo)
             
    -- temperature near the black daisies
    model.tempNearBlack = tempNearDaisy (model.planetTemp, 
            model.planetAlbedo, model.blackAlbedo)
        
    -- indicated growth rate for the white daisies
    model.indWhiteGrowthRate = daisyGrowthRate (model.tempNearWhite)
    
    -- indicated growth rate for the white daisies
    model.indBlackGrowthRate = daisyGrowthRate (model.tempNearBlack)
        
    -- actual growthRate is constained by the empty area of the planet
        
    model.whiteGrowthRate = growthReduction (model.indWhiteGrowthRate, model.emptyArea)
    -- white area                            
    model.whiteArea = model.whiteArea + model.whiteArea * 
                            (model.whiteGrowthRate - model.decayRate) -- + 0.01
        
    -- actual growthRate is constained by the empty area of the planet
        
    model.blackGrowthRate = growthReduction 
                                    (model.indBlackGrowthRate, model.emptyArea)
    -- white area                            
    model.blackArea = model.blackArea + model.blackArea * 
                            (model.blackGrowthRate - model.decayRate) -- + 0.01
        
    model.daisyArea = model.blackArea + model.whiteArea
    model.emptyArea = 1.0 - model.daisyArea
end

DaisyWorld = Model {
    
            -- sun luminosity (this is the main variable of the model)
        -- values beteween 0.70 and 1.6  support life in Daisyworld
        sunLuminosity         =  0.8, 
    
        -- variables 
        -- default values are approximate equilibrium values 
        -- when sun luminosity is equal to 1.0
        -- hint: change the starting conditions when changing the sun luminosity
        -- increase the white area for a strong sun
        -- increase the black area for a faint sun
        planetArea           =   1.0,
        daisyArea            =   0.673,
        whiteArea            =   0.40, 
        blackArea            =   0.273, 
        emptyArea            =   0.327,	
        planetAlbedo         =   0.532,
        planetTemp           =   295.0, -- in K (corresponds to 22 C) 
	    	                                                             
	    -- params (fixed assumptions)
        whiteAlbedo          =   0.75,
        blackAlbedo          =   0.25, 
        soilAlbedo           =   0.50, 
        decayRate            =   0.3,
        
        -- simulation time
        step       =  0.01,
        --updateTime =   100,
        finalTime  =   80,


        init = function (model)
            
            model.timer = Timer {
                    Event { 
                        action = function (event)
                                    growDaisies(model)
                                    --if ((event:getTime() % 100) == 0 ) then 
                                    --    model.sunLuminosity =  model.sunLuminosity + 0.01 
                                    --end
                                    model:notify()

                                end
                        }
            }

        end
}
---------------------------------------------------------------------------------------
-- creates an instance of the model
myDaisyWorld = DaisyWorld { sunLuminosity = 0.7 }

myDaisyWorld.graphics = {
    
    graphics1 = Chart {
                target = myDaisyWorld,
                select = {"blackArea", "whiteArea", "emptyArea"}
    },
    graphics2 = Chart {
                target = myDaisyWorld,
                select = {"planetAlbedo"}
    },
    graphics3 = Chart {
                target = myDaisyWorld,
                select = {"sunLuminosity"}
    },
    graphics4 = Chart {
                target = myDaisyWorld,
                select = {"planetTemp"}
    },
    graphics5 = Chart {
                target = myDaisyWorld,
                select = {"daisyArea"}
    }                        
}
myDaisyWorld:execute ()