import numpy as np
import scipy
import scipy.io
import pylab
import numpy
import glob
import pyfits

def mklc(t, nspot=200, incl=(scipy.pi)*5./12., amp=1., tau=30.5, p=10.0):

    diffrot = 0.

    ''' This is a simplified version of the class-based routines in
    spot_model.py. It generates a light curves for dark, point like
    spots with no limb-darkening.

    Parameters:
    nspot = desired number of spots present on star at any
            one time
    amp = desired light curve amplitude
    tau = characteristic spot life-time
    diffrot = fractional difference between equatorial and polar
              rotation period
    (unit of time is equatorial rotation period)'''

    # print('Period = ', p)
    dur = (max(t) - min(t))

    # (crude estimate of) total number of spots needed during entire
    # time-series
    nspot_tot = int(nspot * dur / 2 / tau)

    # uniform distribution of spot longitudes
    lon = scipy.rand(nspot_tot) * 2 * scipy.pi

   # distribution of spot latitudes uniform in sin(latitude)
    lat = scipy.arcsin(scipy.rand(nspot_tot))

    # spot rotation rate optionally depends on latitude
    period = ((scipy.sin(lat) - 0.5) * diffrot + 1.0 ) * p
    period0 = scipy.ones(nspot_tot) * p

    # all spots have the same maximum area
    # (crude estimate of) filling factor needed per spot
    ff = amp / scipy.sqrt(nspot)
    scale_fac = 1
    amax = scipy.ones(nspot_tot) * ff * scale_fac

    # all spots have the evolution timescale
    decay = scipy.ones(nspot_tot) * tau

    # uniform distribution of spot peak times
    # start well before and end well after time-series limits (to
    # avoid edge effects)
    extra = 3 * decay.max()
    pk = scipy.rand(nspot_tot) * (dur + 2 * extra) - extra

    # COMPUTE THE LIGHT CURVE
    # print("Computing light curve...")
    time = numpy.array(t - min(t))

    area_tot = scipy.zeros_like(time)
    dF_tot = scipy.zeros_like(time)
    dF_tot0 = scipy.zeros_like(time)

    # add up the contributions of individual spots
    for i in range(nspot_tot):

        # Spot area
        if (pk[i] == 0) + (decay[i] == 0):
            area = scipy.ones_like(time) * amax[i]
        else:
            area = amax[i] * \
                scipy.exp(-(time - pk[i])**2 / 2. / decay[i]**2)
        area_tot += area

        # Fore-shortening
        phase = 2 * scipy.pi * time / period[i] + lon[i]
        phase0 = 2 * scipy.pi * time / period0[i] + lon[i]
        mu = scipy.cos(incl) * scipy.sin(lat[i]) + \
            scipy.sin(incl) * scipy.cos(lat[i]) * scipy.cos(phase)
        mu0 = scipy.cos(incl) * scipy.sin(lat[i]) + \
            scipy.sin(incl) * scipy.cos(lat[i]) * scipy.cos(phase0)
        mu[mu < 0] = 0.0
        mu0[mu0 < 0] = 0.0

        # Flux
        dF_tot -= area * mu
        dF_tot0 -= area * mu0

    amp_eff = dF_tot.max()-dF_tot.min()
    nspot_eff = area_tot / scale_fac / ff

    res0 = scipy.array([nspot_eff.mean(), ff, amp_eff])
    res1 = scipy.zeros((4, len(time)))

    res1[0,:] = time
    res1[1,:] = area_tot
    res1[2,:] = dF_tot
    res1[3,:] = dF_tot0

    # print('Used %d spots in total over %d rotation periods.' % (nspot_tot, dur))
    # print('Mean filling factor of individual spots was %.4f.' % ff)
    # print('Desired amplitude was %.4f, actual amplitude was %.4f.' \
    #         % (amp, amp_eff))
    # print('Desired number of spots at any one time was %d.' % nspot)
    return res0, res1

[((1870, 1892), 'scipy.zeros_like', 'scipy.zeros_like', (['time'], {}), '(time)\n', (1886, 1892), False, 'import scipy\n'), ((1906, 1928), 'scipy.zeros_like', 'scipy.zeros_like', (['time'], {}), '(time)\n', (1922, 1928), False, 'import scipy\n'), ((1943, 1965), 'scipy.zeros_like', 'scipy.zeros_like', (['time'], {}), '(time)\n', (1959, 1965), False, 'import scipy\n'), ((1055, 1076), 'scipy.rand', 'scipy.rand', (['nspot_tot'], {}), '(nspot_tot)\n', (1065, 1076), False, 'import scipy\n'), ((1208, 1229), 'scipy.ones', 'scipy.ones', (['nspot_tot'], {}), '(nspot_tot)\n', (1218, 1229), False, 'import scipy\n'), ((1350, 1367), 'scipy.sqrt', 'scipy.sqrt', (['nspot'], {}), '(nspot)\n', (1360, 1367), False, 'import scipy\n'), ((1494, 1515), 'scipy.ones', 'scipy.ones', (['nspot_tot'], {}), '(nspot_tot)\n', (1504, 1515), False, 'import scipy\n'), ((933, 954), 'scipy.rand', 'scipy.rand', (['nspot_tot'], {}), '(nspot_tot)\n', (943, 954), False, 'import scipy\n'), ((1397, 1418), 'scipy.ones', 'scipy.ones', (['nspot_tot'], {}), '(nspot_tot)\n', (1407, 1418), False, 'import scipy\n'), ((1698, 1719), 'scipy.rand', 'scipy.rand', (['nspot_tot'], {}), '(nspot_tot)\n', (1708, 1719), False, 'import scipy\n'), ((2132, 2153), 'scipy.ones_like', 'scipy.ones_like', (['time'], {}), '(time)\n', (2147, 2153), False, 'import scipy\n'), ((2225, 2278), 'scipy.exp', 'scipy.exp', (['(-(time - pk[i]) ** 2 / 2.0 / decay[i] ** 2)'], {}), '(-(time - pk[i]) ** 2 / 2.0 / decay[i] ** 2)\n', (2234, 2278), False, 'import scipy\n'), ((2455, 2470), 'scipy.cos', 'scipy.cos', (['incl'], {}), '(incl)\n', (2464, 2470), False, 'import scipy\n'), ((2473, 2490), 'scipy.sin', 'scipy.sin', (['lat[i]'], {}), '(lat[i])\n', (2482, 2490), False, 'import scipy\n'), ((2545, 2561), 'scipy.cos', 'scipy.cos', (['phase'], {}), '(phase)\n', (2554, 2561), False, 'import scipy\n'), ((2576, 2591), 'scipy.cos', 'scipy.cos', (['incl'], {}), '(incl)\n', (2585, 2591), False, 'import scipy\n'), ((2594, 2611), 'scipy.sin', 'scipy.sin', (['lat[i]'], {}), '(lat[i])\n', (2603, 2611), False, 'import scipy\n'), ((2666, 2683), 'scipy.cos', 'scipy.cos', (['phase0'], {}), '(phase0)\n', (2675, 2683), False, 'import scipy\n'), ((1150, 1164), 'scipy.sin', 'scipy.sin', (['lat'], {}), '(lat)\n', (1159, 1164), False, 'import scipy\n'), ((2507, 2522), 'scipy.sin', 'scipy.sin', (['incl'], {}), '(incl)\n', (2516, 2522), False, 'import scipy\n'), ((2525, 2542), 'scipy.cos', 'scipy.cos', (['lat[i]'], {}), '(lat[i])\n', (2534, 2542), False, 'import scipy\n'), ((2628, 2643), 'scipy.sin', 'scipy.sin', (['incl'], {}), '(incl)\n', (2637, 2643), False, 'import scipy\n'), ((2646, 2663), 'scipy.cos', 'scipy.cos', (['lat[i]'], {}), '(lat[i])\n', (2655, 2663), False, 'import scipy\n')]