%PROGRAM    readSEM2
%PURPOSE    The purpose of this program is to read in the archive SEM2
%           files provided by SEC and available on the web at NGDC (full resoultuion files). 
%           It reads in the archive data using the pogram called
%           unpack_sem2 and passes back data,columns and units.
%           unpack_sem2 is based on the unpack program provided by
%           Dave Evans at SEC. unpack_sem2.f simply includes the mex file
%           wrapper so that the unpack routine can be clled from matlab.
%USAGE      [data,columns,units]=readSEM2_new('filename')
%           where data is the returned data. The first column is the time
%           given in matlab datenumber. All the columns are defined by the
%           cell array, columns.
%DATE       10/2004
%Notes      WARNING: RIGHT NOW THE PROGRAM ONLY OUTPUTS MEPED and ORBIT DATA!!!
%           I HAVE NOT DEALT WITH THE TED DATA WHICH IS ALSO INCLUDED IN THE RAW
%           DATA FILES
%           The MEPED SEM-2 system measures electrons in 3 integral
%           channels in two look directions.
%           e1  >30keV  telescope detector  1 sec accumulation
%           e2  >100keV telescope detector  1 sec accumulation
%           e3  >300keV telescope detector  1 sec accumulation  
%           The electron measurements are sensitive to protons between 135 keV and 1MeV.
%           The system measures protons in 5 differential energy bands and 3 integral dome
%           detector channels;
%           p1  30-80keV        telescope detector  1 sec accumulation
%           p2  80-240keV       telescope detector  1 sec accumulation
%           p3  240-800keV      telescope detector  1 sec accumulation
%           p4  800-2500keV,    telescope detector  1 sec accumulation
%           p5  2.5-6.9 MeV     telescope detector  16 sec accumulation
%           p6  >6.9 MeV        telescope detector
%           p6  >16 MeV         dome detector        2 sec accumulation
%           p7  >36             dome detector        2 sec accumulation
%           p8  >80 MeV.        dome detector        2 sec accumulation
%           The proton telescope detectors are sensitive to >1MeV electrons
%           Even though the protons have differntial
%           measurements the units are still #/cm^2-s-str because the dE's are not
%           given. I've kept the dome channels in units of counts/seconds
%           because there's some ambiguity in translating them to flux
%           measurements. The returned time is given at 1 second resolution but all
%           energies and look direction are not measured every second.
%           First the 0 degree particles are measured and accumulated for
%           1 second and then the 90 degree particles are measured and
%           accumulated for 1 second. The ominidirectional protons are
%           measured once every every 2 seconds and accumulated for 2
%           seconds. All times are tagged at the end of the accumulation
%           period!



function [data,columns,units]=readSEM2_new(infile)

[time1,mep1,orb1,stat1,count1]=unpack_sem2_safe(infile);
count=double(count1)-1;
ttime=double(time1(1:count,:));

%calculate matlab time from date vector
dtest=datenum((ttime(:,1)),0,(ttime(:,2)),0,0,0);
[TY,MO,D,S,S,S]=datevec(dtest);
%[MO,D,S,S,S]=date((ttime(:,1)),(ttime(:,2)));
tempmatlabd=datenum((ttime(:,1)),MO,D,(ttime(:,3)),(ttime(:,4)),(ttime(:,5)));
%At the end of the files, time is sometimes funny. Often there is data from
%the previous day or repeat records. This gets rid of the funny data
%+++++++++++++++++Correct the time++++++++++++
goodrec1=find(tempmatlabd>=tempmatlabd(1));
[temp,I]=sort(tempmatlabd(goodrec1));
temp2=diff(temp);
ind1=find(temp2~=0);
ind=I(ind1);
goodrec=goodrec1(I(ind1));
matlabd=tempmatlabd(goodrec);
time=ttime(goodrec,:);
clear ttime tempmatlabd

mep=mep1(goodrec,:);
clear mep1 ;
mep(mep<0)=NaN;

orb=orb1(goodrec,:);
clear orb1 ;

stat=stat1(goodrec,:);
clear stat1 ;

data=zeros(length(time),39).*NaN;

%I'm setting this up so that the columns are similar to SEM1
data(:,1)=matlabd;

data(:,2:10)=mep(:,10:18); %This is the field aligned component
data(:,11:19)=mep(:,1:9); %This is perpindicular component
data(:,20:23)=mep(:,19:22); %This is omnidirectional data

data(:,24)=stat(:,3);%This is meped inflight cal flag
data(:,25)=orb(:,2);%lat
data(:,26)=orb(:,3);%lon
data(:,27)=time(:,6);%alti
data(:,28)=time(:,7);%orbit
data(:,29:32)=orb(:,4:7);%Br,Be,Bs,Btot local
data(:,33:34)=orb(:,10:11);%0,90 pa local
data(:,35:36)=orb(:,18:19);%0,90 pa local
data(:,37)=orb(:,24);%L
data(:,38)=orb(:,25);%maglat
data(:,39)=orb(:,26);%MLT
%=============Now do the calibration of the MEPED data
data(:,2:19)=data(:,2:19).*100; %This puts all the directional detectors in #/cm2-s-str
%For the domes, the goemetric factor depends on the pitch angle so I left
%it as counts. The only thing to consider is the integration time
data(:,20:21)=data(:,20:21)./2; 
data(:,22:23)=data(:,22:23)./4; 

columns={'time',
  'proton-1 field aligned .03-.08 MeV',
  'proton-2 field aligned .08-.25 MeV',
   'proton-3 field aligned .25-.80 MeV',
   'proton-4 field aligned .80-2.50 MeV',
   'proton-5 field aligned 2.50-6.90 MeV',
   'proton-6 field aligned >6.00 MeV',
'electron-1 field aligned .03-2.50 MeV',
'electron-2 field aligned .10-2.50 MeV',
'electron-3 field aligned .30-2.50 MeV',
'proton-1 perpendicular .03-.80 MeV' ,
'proton-2 perpendicular .80-.24 MeV',
'proton-3 perpendicular .24-.80 MeV',
'proton-4 perpendicular .80-2.50 MeV',
'proton-5 perpendicular 2.50-6.90 MeV',
'proton-6 perpendicular aligned >6.00 MeV',
'electron-1 perpendicular .03-2.50 MeV',
'electron-2 perpendicular .10-2.50 MeV',
'electron-3 perpendicular .30-2.50 MeV',
'proton-6 dome 16 MeV',
'proton-7 dome 35 MeV',
'proton-8 dome 80 MeV',
'proton-9 dome 140 MeV',
'In-flight Calibration',
'lat (-90 to 90)',
 'lon (0-360)',
 'altitude',
 'orbit',
 'Br IGRF',
 'Be IGRF',
 'Bs IGRF',
 'Btot IGRF',
 'pitch angle (local) field aligned detector',
 'pitch angle (local) perpendicualar detector',
 'pitch angle (foot of the field) field aligned detector',
 'pitch angle (foot of the field) perpendicular detector',
 'L IGRF',
 'magnetic latitude',
 'LT (0-360)'};

units={'matlab datenum',
  '#/cm^2-s-str',
  '#/cm^2-s-str',
   '#/cm^2-s-str',
   '#/cm^2-s-str',
   '#/cm^2-s-str',
   '#/cm^2-s-str',
'#/cm^2-s-str',
'#/cm^2-s-str',
'#/cm^2-s-str',
'#/cm^2-s-str' ,
'#/cm^2-s-str',
'#/cm^2-s-str',
'#/cm^2-s-str',
'#/cm^2-s-str',
'#/cm^2-s-str',
'#/cm^2-s-str',
'#/cm^2-s-str',
'#/cm^2-s-str',
'#/cm^2-s-str',
'#/cm^2-s',
'#/cm^2-s',
'#/cm^2-s',
' ',
'degrees -90-90',
 'degrees(0-360)',
 'km',
 '#',
 'nT',
 'nT',
 'nT',
 'nT',
 'degrees',
 'degrees',
 'degrees',
 'degrees',
 'L IGRF',
 'degrees',
 'degrees (0-360)'};