Generates temperature maps with manually chosen temperatures, sizes, and positions of star spots then writes the artificially observed flux and phase data into a file. Written for Julia 0.6.4. Repository ROTIR.jl by Fabien Baron has been updated for Julia 1.0.1 and will no longer be compatible with this script.

ControlData.jl

#Generate fake light-curve
include("lci.jl");
include("lciplot.jl");
include("geometry.jl");
include("oiplot.jl");
include("oistars.jl");
using PyPlot
using PyCall
using ROTIR.jl

#fake star params
nepochs = 400;
phase = linspace(0,20,nepochs);

stellar_parameters = Array{starparameters}(nepochs);
params = [1.,    # radius (fixed for LCI)
          5200., # temperature
          0.,    # fractional critical velocity (assume 0 for non-rapid rotators)
          60.0,  # inclination angle
          0.0,   # position angle
          0,     # rotation angle (for LCI we use phase instead)
          3.0,   # rotation period (for Kepler LCI, does not matter)
          [3,0.2],   # Hestroffer limb darkening coefficient
          0.08,  # von Zeipel exponent
          0,     # differential rotation coefficient 1
          0      # differential rotation coefficient 2
          ];

#period = 3
for i=1:nepochs
    stellar_parameters[i]=starparameters(params[1],params[4],params[5],360.0*(1/3)*phase[i],params[7],params[8],params[10],params[11]);
end

n = 4; #Healpix tesselation level
star_epoch_geom = create_geometry(healpix_round_star(n), stellar_parameters);
npix = star_epoch_geom[1].npix;
polyflux = setup_lci(star_epoch_geom);

### W A R N I N G: will generate [nepochs/(2*maximum(phase))] images
include("ControlMap.jl");
#returns x_all

flux=[]
fluxerr=[]
fl=zeros(nepochs);
for i=1:length(nmaps)
    fl=polyflux*x_all[i,:]
    flerr = 5e-4*maximum(fl)*ones(nepochs)
    fl += flerr.*randn(nepochs);
    for j=fl
        flux=push!(flux, j)
    end
    for k=flerr
        fluxerr=push!(fluxerr, k)
    end
end
lcidata  = LCI([], phase, flux, fluxerr, nepochs);
write_lci("FAKE_flux.txt", lcidata)

ControlMap.jl

using ROTIR.jl
using LibHealpix

neighbors,south_neighbors,west_neighbors,south_neighbors_reverse,west_neighbors_reverse = tv_neighbours_healpix(n);

function circular_spots(x_phot::Array{Float64,1}, ipix1::Int64, temp1::Float64, ipix2::Int64, temp2::Float64)

    h=ipix->vcat(neighbors[ipix]...)
        function neighlvl(ipix1, n)
            list = h(ipix1);
            for k=1:n-1
            list = h(list)
            end
         return list
        end

        x = deepcopy(x_phot);

            x[ipix1] = temp1
            x[neighlvl(ipix1,2)] = temp1
            y = neighlvl(ipix1,2)
                for u = 1:length(y)
                ipixN = south_neighbors_reverse[y[u]]
                ipixS = south_neighbors[y[u]]
                x[ipixN] = temp1
                x[ipixS] = temp1
                x[neighbors[south_neighbors[south_neighbors[ipix1]]][8]] = 5200.
                x[neighbors[south_neighbors[south_neighbors[ipix1]]][2]] = 5200.
                x[neighbors[south_neighbors[south_neighbors[ipix1]]][1]] = 5200.
                x[neighbors[neighbors[neighbors[ipix1][4]][5]][5]] = 5200.
                x[neighbors[neighbors[neighbors[ipix1][5]][5]][6]] = 5200.
                x[neighbors[neighbors[neighbors[ipix1][5]][5]][5]] = 5200.
            end

            x[ipix2] = temp2
            x[neighlvl(ipix2,1)] = temp2
            y = neighlvl(ipix2,1)
                for v = 1:length(y)
                ipixN = south_neighbors_reverse[y[v]]
                ipixS = south_neighbors[y[v]]
                x[ipixN] = temp2
                x[ipixS] = temp2
                x[south_neighbors[south_neighbors[ipix2]]] = 5200.
                x[south_neighbors[south_neighbors[ipix2]]] = 5200.
                x[south_neighbors[south_neighbors[ipix2]]] = 5200.
                x[neighbors[neighbors[ipix2][4]][5]]= 5200.
                x[neighbors[neighbors[ipix2][5]][5]] = 5200.
                x[neighbors[neighbors[ipix2][5]][5]] = 5200.
                end
            return x
end

nmaps = Int64(nepochs/(2*maximum(phase))) #change of spots per day (2x)
x_phot = 5200*ones(npix);
x_all = ones(nmaps, length(x_phot)); # collects temp. maps at all epochs
ipix1 = 300; #center pixels of spots:
ipix2 = 2300;
for i = 1:nmaps
    ipix1 = west_neighbors[ipix1]; #motion spot 1
    ipix1 =Int64(round(ipix1))
    ipix2 = west_neighbors[ipix2]; #motion spot 2
    ipix2 =Int64(round(ipix2))
    x_all[i,:] =  circular_spots(x_phot, ipix1, 4200., ipix2, 4500.);
end

#plot temperature maps
for l=1:nmaps
       mollplot_temperature(x_all[l,:]); PyPlot.draw();PyPlot.pause(1.0);
end