function [LC,ZP,BinLC,BinConLC]=read_mark_lc(File,BinVec,BinType,BinSize,Randomize,RemoveDuplicateLines)
% Read PTF subtraction light curves generated by Mark S. program.
% Package: VO.PTF
% Description: Read PTF subtraction light curves generated by Mark S. program.
% Input  : - File.
%          - Binning vector [MJD]. This is a column vector of edges defines
%            the bins.
%            Alternatively, a two column vector in which the first (second)
%            column is for the start (end) time of each bin.
%          - In addition to a predefined bins, the program can calculate
%            the light curve in a constant bins of either constant
%            Number of points or constant time span.
%            BinType can be:
%            'Time' - for constant time span.
%            'Num'  - for constant number of points.
%          - Time span or number of points in bins.
%          - Randomize the data points (i.e., assign random dates),
%            {true|false}. Default is false.
%          - A flag indicating if to remove lines with the same exact JD.
%            Default is true.
% Output : - Light curve [MJD, Mag, MagErr, 3sigma limit, Flux, Flux error]
%          - Zero point
%          - Binned light curve, in the pre-defined bins defined
%            in BinVec. The matrix contains the following columns:
%            (1) Mean time [MJD]
%            (2) Err Time low [day]
%            (3) Err Time up [day]
%            (4) Number of observations
%            (5) Mag
%            (6) Mag Err
%            (7) Mag limit
%            (8) median Flux
%            (9) Flux Err
%            (10) Mag Err derived by bootstrap simulations
%            (11) Flux Err derived by bootstrap simulations
%          - Binned light curve in the constant number/time-span bins.
%            The matrix contains the following columns:
%            (1) Mid time
%            (2) Median flux in bin.
%            (3) Std of flux in bin divided by the number of points in bin.
%            (4) Median flux error in bin.
%            (5) Number of measurements in bin.
%            (6) Bootstrap Std of flux in bin.
%            (7) Jackknife Std of flux in bin.
% Tested : Matlab 7.13
%     By : Eran O. Ofek / Aviram Steinbok  Apr 2012
%    URL : http://weizmann.ac.il/home/eofek/matlab/
%-------------------------------------------------------------------------------

Def.BinVec    = [];
Def.BinType   = [];
Def.BinSize   = [];
Def.Randomize = false;
Def.RemoveDuplicateLines = true;
if (nargin==1)
    BinVec         = Def.BinVec;
    BinType        = Def.BinType;
    BinSize        = Def.BinSize;
    Randomize = Def.Randomize;
    RemoveDuplicateLines = Def.RemoveDuplicateLines;
elseif (nargin==2)
    BinType        = Def.BinType;
    BinSize        = Def.BinSize;
    Randomize = Def.Randomize;
    RemoveDuplicateLines = Def.RemoveDuplicateLines;
elseif (nargin==3)
    BinSize        = Def.BinSize;
    Randomize = Def.Randomize;
    RemoveDuplicateLines = Def.RemoveDuplicateLines;
elseif (nargin==4)
    Randomize = Def.Randomize;
    RemoveDuplicateLines = Def.RemoveDuplicateLines;
elseif (nargin==5)
    RemoveDuplicateLines = Def.RemoveDuplicateLines;
else
    % do nothing
end

%JDoffset = -2400000.5 +2450000;
JDoffset = 0;

BootStrapFun = 'mean';
B = 1e-10;  % lapyitude sofetning paramerter

Nsim = 100;

Nsig = 3;  % sigma for background limit

FID = fopen(File,'r');
line=Util.IO.file2str(FID,'cell');
fclose(FID);
hdrlines=find(Util.cell.isempty_cell(strfind(line,'MJD counts'))==0);
if isempty(hdrlines)
    hdrlines=12; % use 12 as default 
end
FID = fopen(File,'r');
C = textscan(FID,'%f %f %f %f %*[^\n]','Headerlines',hdrlines,'CommentStyle','#');
fclose(FID);

if (RemoveDuplicateLines)
    [C{1},ind]=unique(C{1});
    C{2}=C{2}(ind);
    C{3}=C{3}(ind);
    C{4}=C{4}(ind);
end

%if (isempty(C{1})),
%    FID = fopen(File,'r');
%    C = textscan(FID,'%f %f %f %f %*[^\n]','Headerlines',7,'CommentStyle','#');
%    fclose(FID);
%end

Lup = AstroUtil.spec.luptitude(C{2},1,B);
%LC = [C{1}, C{4}-2.5.*log10(C{2}), 1.086.*C{3}./C{2}, C{4}-2.5.*log10(Nsig.*C{3}), C{2}, C{3}];

LC = [C{1}, C{4}+Lup, 1.086.*C{3}./C{2}, C{4}-2.5.*log10(Nsig.*C{3}), C{2}, C{3}];

LC(:,1) = LC(:,1) - JDoffset;
if (Randomize)
   Nd = size(LC,1);
   RandInd = randperm(Nd,Nd);
   LC = [LC(:,1),LC(RandInd,2:end)];
end

ZP = C{4};
BinLC    = [];
BinConLC = [];

if (~isempty(BinVec))
   if (size(BinVec,2)==2)
      % two column matrix [start end] of each bin
      Nbin = size(BinVec,1);
      Ncol = 2;
   else
      Nbin  = length(BinVec) - 1;
      Ncol  = 1;
   end
   BinLC = zeros(Nbin,11);
   for Ibin=1:1:Nbin
     switch Ncol
       case 1
          Isel   = find(LC(:,1)>=BinVec(Ibin) & LC(:,1)<BinVec(Ibin+1));
       case 2
          Isel   = find(LC(:,1)>=BinVec(Ibin,1) & LC(:,1)<BinVec(Ibin,2));
       otherwise
          error('Illegal Ncol option');      
      end
      N      = length(Isel);
      MeanT  = mean(LC(Isel,1));
      if (N==0)
         MinT = NaN;
         MaxT = NaN;
         ErrTl= NaN;
         ErrTu= NaN;
      else
         MinT   = min(LC(Isel,1));
         MaxT   = max(LC(Isel,1));
         ErrTl  = MeanT - MinT;
         ErrTu  = MaxT  - MeanT;
      end
      [Flux,FluxErr] = Util.stat.wmean(LC(Isel,5:6));
      
      Flux = median(LC(Isel,5));
      
      MeanZP = mean(ZP(Isel));
      Mag    = MeanZP - 2.5.*log10(Flux);
      % Mag = convert.luptitude(Flux,10.^(0.4.*MeanZP));
      Err    = 1.086.*FluxErr./Flux;
      MagLimit = MeanZP - 2.5.*log10(FluxErr.*Nsig);

      % bootstrap for errors
      % select N out of N with repeations
      Nsel  = length(Isel);
      if (Nsel>1)
          MagStD  = Util.stat.bootstrap_std(LC(Isel,2),BootStrapFun,Nsim.*Nsel);
          FluxStD = Util.stat.bootstrap_std(LC(Isel,5),BootStrapFun,Nsim.*Nsel);
      else
          MagStD  = NaN;
          FluxStD = NaN;
      end
      %Nsel  = length(Isel);
      %ThetaMag = zeros(Nsim,1);
      %ThetaErr = zeros(Nsim,1);
      %ThetaMagLimit = zeros(Nsim,1);
      %for Isim=1:1:Nsim,
      %   Irand = ceil(rand(Nsel,1).*Nsel);
      %   [ThetaFlux,ThetaFluxErr] = wmean(LC(Isel(Irand),5:6));
      %   ThetaMeanZP = mean(ZP(Isel(Irand)));
      %   ThetaMag(Isim)      = ThetaMeanZP - 2.5.*log10(ThetaFlux);
      %   ThetaErr(Isim)      = 1.086.*ThetaFluxErr./ThetaFlux;
      %   ThetaMagLimit(Isim) = ThetaMeanZP - 2.5.*log10(ThetaFluxErr.*Nsig);
      %end
      %MeanThetaMag = mean(ThetaMag);
      %MeanThetaErr = mean(ThetaErr);
      %MeanThetaMagLimit = mean(ThetaMagLimit);
      %MagStD = sqrt(sum((ThetaMag - MeanThetaMag).^2)./(Nsim-1));
      %--- end of bootstrap section ---

      BinLC(Ibin,:) = [MeanT,ErrTl, ErrTu,N,Mag,Err,MagLimit,Flux,FluxErr, MagStD, FluxStD];
   end

   % remove entries with 0 observations
   %BinLC = BinLC(BinLC(:,4)>0,:);
end


if (nargin>2 && nargout>3)
    % calculate LC in constant bins
    switch lower(BinType)
        case 'time'
            N = size(LC,1);
            BinConLC = zeros(N,7);
            for I=1:1:N
                Ibin    = find(LC(:,1)>=LC(I,1) & LC(:,1)<(LC(I,1)+BinSize));
                MidTime = 0.5.*(max(LC(Ibin,1)) + min(LC(Ibin,1)));
                MedFlux = median(LC(Ibin,5));
                Npoint  = numel(Ibin);
                StdFlux = std(LC(Ibin,5))./sqrt(Npoint);
                MedErr  = median(LC(Ibin,6));
                bootstrapStdFlux = bootstrap_std(LC(Ibin,5),BootStrapFun,Nsim.*Npoint);
                jackknifeStdFlux = jackknife(LC(Ibin,5),BootStrapFun);
                BinConLC(I,:) = [MidTime, MedFlux, StdFlux, MedErr, Npoint, bootstrapStdFlux, jackknifeStdFlux];

            end

                
        case 'num'
            MidTime = LC(:,1);
            MedFlux = medfilt1(LC(:,5),BinSize);
            StdFlux = stdfilt1(LC(:,5),BinSize)./sqrt(BinSize);
            MedErr  = medfilt1(LC(:,6),BinSize);
            bootstrapStdFlux=zeros(size(LC(:,5)));
            jackknifeStdFlux=zeros(size(LC(:,5)));
            for I=1:1:size(LC(:,5))
                if (((I-(ceil(BinSize/2)))>=0))
                    if I > ((size(LC(:,5))-(ceil(BinSize/2)-ceil(ceil(BinSize/2)-(BinSize/2)))))
                        bootstrapStdFlux(I) = bootstrap_std(LC((I-ceil(BinSize/2)+1):size(LC(:,5)),5),BootStrapFun,Nsim.*BinSize);
                        jackknifeStdFlux(I) = jackknife(LC((I-ceil(BinSize/2)+1):size(LC(:,5)),5),BootStrapFun);                     
                    else
                        bootstrapStdFlux(I) = bootstrap_std(LC((I-ceil(BinSize/2)+1):(I+ceil(BinSize/2)-ceil(ceil(BinSize/2)-BinSize/2)),5),BootStrapFun,Nsim.*BinSize);
                        jackknifeStdFlux(I) = jackknife(LC((I-ceil(BinSize/2)+1):(I+ceil(BinSize/2)-ceil(ceil(BinSize/2)-BinSize/2)),5),BootStrapFun);
                    end
                else
                    bootstrapStdFlux(I) = bootstrap_std(LC(1:(I+ceil(BinSize/2)-ceil(ceil(BinSize/2)-BinSize/2)),5),BootStrapFun,Nsim.*BinSize);
                    jackknifeStdFlux(I) = jackknife(LC(1:(I+ceil(BinSize/2)-ceil(ceil(BinSize/2)-BinSize/2)),5),BootStrapFun);
                end
            end
            Npoint  = BinSize.*ones(size(MidTime));
            BinConLC = [MidTime, MedFlux, StdFlux, MedErr, Npoint, bootstrapStdFlux, jackknifeStdFlux];
        otherwise
            error('Unknwon BinType option');
    end
end

