from __future__ import print_function, division

import sys
import numpy as np
import numpy.random as rand
import random
import matplotlib.pyplot as plt
from pymacula import Star, Spot, MaculaModel
#from ._macula import maculamod as macula


def simulate_LC(Prot, incl, nspots):
  dt = 48.
  Tmax = 1200.
  ndata = 57600.                  # 1200 days of continuous photometry, with resolution of 30 minutes => 1200*24*2 = 57600
  # array of timesteps at which data is "collected"
  ts = np.linspace(dt, Tmax, ndata)

  # parameter vector for the star's intrinsic parameters: inclination, rotation period etc.
  star = Star(incl=incl,  Peq=Prot, kappa2=0.2, kappa4=0.2, 
              c1=0.3999, c2=0.4269, c3=-0.0227,  c4=-0.0839,
              d1=0.3999, d2=0.4269, d3=-0.0227,  d4=-0.0839)

  spots = []                                # list of all spots
  theta_spot = np.zeros((8, nspots))        # parameters of each spot

  mean_lifetime = 100.            # mean of the gaussian distribution which gives the lifetimes of spots
  std_lifetime = 30.              # standard deviation of the gaussian distribution which gives the lifetimes of spots

  # parameters of the kth starspot  
  for k in range(nspots):
    # start time of the k-th spot, drawn uniformly at random from the interval [Prot, Tmax]
    Tstart = rand.random()*Tmax
    # lifetime of the k-th spot, randomly drawn from a gaussian distribution
    lifetime = random.gauss(mean_lifetime, std)
    # end time of the k-th spot
    Tend = Tstart + lifetime

    theta_spot[0,k] = rand.random()*np.pi                     # longitude of the k-th spot at time tref(k)
    theta_spot[1,k] = rand.random()*np.pi/2.                  # latitude of the k-th spot at time tref(k)
    theta_spot[2,k] = rand.random()*10.                       # angular spot size of the k-th spot at time tmax(k)
    theta_spot[3,k] = rand.random()*0.5                       # spot to star flux contrast of the k-th spot

    theta_spot[4,k] = rand.random()*lifetime + Tstart         # time at which the k-th spot is largest
    theta_spot[5,k] = lifetime                                # lifetime of the k-th spot in days 
    theta_spot[6,k] = rand.random()*lifetime                  # ingress duration of the k-th spot in days
    theta_spot[7,k] = rand.random()*lifetime                  # egress duration of the k-th spot in days

    # parameter vector for the k-th starspot
    spot = Spot(lon=theta_spot[0,k], lat=theta_spot[1,k],
                alpha_max=theta_spot[2,k], contrast=theta_spot[3,k],
                tmax=theta_spot[4,k], lifetime=theta_spot[5,k],
                ingress=theta_spot[6,k], egress=theta_spot[7,k])
                  
    # add spot to the list of spots
    spots.append(spot)


  # generate the light curve 
  # ??? inst_offsets and blend_factors np.array([1.,1.]) ???
  model = MaculaModel(t=ts, star=star, spots=spots, nspots=nspots)
  return ts, model(ts) 


def main(argv):
  star_Prot = 30.                   # rotation period of the star
  star_incl = 5.                    # inclination of the star
  star_nspots = 60                  # total number of simulated spots
  
  time, flux = simulate_LC(star_Prot, star_incl, star_nspots)

  # add gaussian noise to the flux
  for i in range(len(flux)):
    flux[i] = flux[i] + random.gauss(0.0, 3e-4) 
  
  # plot the light curve
  plt.plot(time, flux)
  plt.xlabel('Time')
  plt.ylabel('Flux')
  plt.show()  

if __name__ == "__main__":
   main(sys.argv[1:])
