function PlotEvsLMA
%NOTICE: First time you use this -- set loaddata=1 and create global
%variables E, LABEL and ZULUTIME at the command line.
%Code for interpreting Sonde sounding data
    %R.Sonnenfeld
    %Version 3.0 -- August 2006
    %Use for interpreting the CG flash at 2020 on Aug 18, 2004
    % Designed to allow analyzing other flashes by changing "minute" variable and
    % putting different "if" statements.

global E
global LABEL; global ZULUTIME
%
% Definitions for E, B and Gxcal, Gzcal etc. 
% are per G.P.Lu's definitions paper Version 2.0 dated 3/2005
% All other physical quantities are in conventional SI units unless otherwise specified
    
%%%% Set options for processing
load_data=0;    %Set to 1 to load_data, set to 0 if data has already been read in.
dedroop=1;      %Set to zero not to dedroop, Set to 1 to dedroop

%%%% Close figures
close all;

%%%% Setup printing options, fonts, and colors for figures in Windows or Linux
draft=1;        %Set to 1 for draft, 0 for final.  There is no time stamp on final figures.
color=1; %Plot in color
%color=0; %Plot in B/W
if color==1
    xc='b-'  ;xtc='b';  %xcolor for line / text
    yc='r-'  ;ytc='r';  %ycolor for line / text
    zc='k-'  ;ztc='k';  %zcolor
    end
if color==0
    xc='k--'  ;  %xcolor
    yc='k-.'  ;  %ycolor
    zc='k-'  ;  %zcolor
    xc=zc; yc=zc;
    xtc='k';ytc='k';ztc='k';
    end

ostype='windows';
if isunix==1 ostype='linux__'; end
disp(strcat('Detected_',ostype,' OS'));

if ostype=='windows'
    fontn='Times';
else
 %   fontn='Utopia';
    fontn='Century Schoolbook L';
 %   fontn='SansSerif';
    end

    
%%%% Specify the data set (time) and file location.
Sondename='Sonde3';
hour=20   ;
minute=20 ;
ZULUTIME=strcat(int2str(hour),int2str(minute));
dir='/home/sprite1_backup/DATA/20040818/MAY05_ASCII_DATA/'; %For data on feynman
%dir='/DATA/20040818/MAY05_ASCII_DATA/'; %For data on sprite1
%dir='./';
datafile=strcat(dir,'Esonde03_08182004_',ZULUTIME,'.csv') 
if load_data==1
    E=load(datafile) ;
end

%%%% Specify plot limits and title
if minute==20 
    secs='20Z'; 
    ssec=20.9; esec=22.03;
    Emin=-inf;          Emax=inf;
    LABEL='LMAvsE20040818'                ;
    titlelabel=strcat(LABEL,ZULUTIME,secs);
    titletext=['LMA vs. E-field',' ',titlelabel];
end

%================================================================
%Location of the balloon at time of analyzed flash, from Esonde03_Location.dat

if minute==20
%**** For 2020 flash
        lat_son=33.864212;  % For the 2020 flash
        lon_son=-107.09716; 
        alt_son=2992.8;   
    %Ground level beneath the balloon, in m (MSL)
        alt_gnd=1813      ;
    %Ground level beneath a particular flash, in m (MSL)
        alt_flash=1768    ; 
    end

%================================================================
%%%% Specify instrument dependent variables
%time constant of the instrument, for dedrooping
    if Sondename=='Sonde3' tau=1.0;  end
    if Sondename=='Sonde5' tau=0.33;  end
    if Sondename=='Sonde6' tau=0.33;  end
    if Sondename=='Sonde7' tau=0.33;  end

%Correction factor, best up to now. 
G_xcal=1.8;G_ycal=1.8;G_zcal=2.5;  %G defined as per GP Lu paper of 3/2005 -- old values
G_xcal=2.7;G_ycal=2.7;G_zcal=2.4;  %Based on reanalysis 8/2006

%Overall charge amplifier gain;
A=0.01;  %Area of plate in square meters
if Sondename=='Sonde3' C=22e-9; end %Front end capacitor, farads 
beta=10; %2nd stage gain [V/V] (dimensionless)
eps0=8.86e-12; %epsilon-0 (farad/meter)
gamma=C/(eps0*beta*A)
gamma=gamma/1000; %Puts it in kV/meter instead of Volts/meter
      
%local declination, from NGDC at NOAA
%http://www.ngdc.noaa.gov/seg/geomag/jsp/struts/calcDeclination
    declination=10.05 ;
    inclination=61.395;

 
%Sample E-field data file
%Epoch      Microsec    V0         V1        V2        V3        PPS       Bx          By       Bz     Error
%1092859959, 39705884, 0.010986, 0.156250, 0.064087, 0.344543, 0.021362, 0.418701, -1.052246, 1.746521, 00
%1092859959, 39705984, 0.010681, 0.155029, 0.064087, 0.344543, 0.023193, 0.418091, -1.050568, 1.747131, 00
%1092859959, 39706084, 0.010681, 0.156250, 0.065308, 0.346680, 0.023193, 0.416361, -1.051127, 1.747335, 00
%1092859959, 39706284, 0.010681, 0.155945, 0.063171, 0.345154, 0.020142, 0.417175, -1.050293, 1.747131, 00

%Set up the time scale for Sonde data
        time_son_all=E( : ,2)/1E6;
        len_son=length(E);
        s_index=find(time_son_all>=ssec,1); %Find first Sonde data point that exceeds start-time ssec
        %This only enables processing to go faster, because end up working
        %with smaller data arrays.
        
        %Only use data points that exceed start-time
        time_son=E(s_index:len_son,2)/1E6;
        off0=-0.0031;  %These offsets were obtained by taking the median of 1 minute
        off1=0.1492; %of data from the 8/18/2004 2021 data set
        off2=0.1056; %They are specific to the 8/18 flight
        off3=0.3436; %Comparing to 2049 data set, all the offsets drifted down by approx 30 mV.
    
%Electric field data
        V0=E(s_index:len_son ,3)-off0;  
        V1=E(s_index:len_son ,4)-off1;
        V2=E(s_index:len_son ,5)-off2;
        V3=E(s_index:len_son ,6)-off3; 
        
        PPS=E(s_index:len_son ,7);
%Magnetic Field data
        SBx=E(s_index :len_son ,8)/4.0;  %B-field data was already smooothed!!
        SBy=E(s_index :len_son ,9)/4.0;
        SBz=E(s_index :len_son ,10)/4.0;
        
        SBx_all=E( : ,8)/4.0;  %B-field data was already smooothed!!
        SBy_all=E( : ,9)/4.0;
        SBz_all=E( : ,10)/4.0;
        
        SBz_filt=lowpass(SBz_all,10000,0.1); 
        SBz_max=max(SBz_filt)
        k=cot(pi*inclination/180);
        
        dtheta_B=(180/pi)*k*(SBz-SBz_max)./(SBz_max);
        if Sondename == 'Sonde3'
            SBy=-1.0*SBy; SBy_all=-1.0*SBy_all;
            SBz=-1.0*SBz; SBz_all=-1.0*SBz_all;
        end  %This fixes the backward B-sensors problem on 2004 data
        dtheta_B_all=(180/pi)*k*(SBz_all-SBz_max)./(SBz_max);
        
% Calculate the clockwise angle between magnetic field and channel 0.
        phi_B=atan2d(SBy,SBx);  %phi_B is in degrees
        phi_B_all=atan2d(SBy_all,SBx_all);  
%        d_angle=mod(phi_B+180,360);
%        s=lowpass(phi_B,10000,0.02);  That's about the best compromise
%        between attenuating the 30 Hz noise and passing an accurate
%        representation of the data.  

        dt=1E-4;  %Data rate is 10000  points/second

        V0_d=V0; V1_d=V1; V2_d=V2; V3_d=V3;
   if dedroop==1
            V0_d=V0+dt*cumsum(V0)/tau;  % The _d suffix on the variable means dedrooped
            V1_d=V1+dt*cumsum(V1)/tau;
            V2_d=V2+dt*cumsum(V2)/tau;
            V3_d=V3+dt*cumsum(V3)/tau;
    %Note that dedrooping doesn't really make sense with gaps in data.
    %dedrooping should be restarted at each new data packet.  At the moment
    %we aren't doing this.  Thus the overall offset of each new data packet
    %has no physical reality.  However the change of each packet works
    %fine.
    else
        warning('Informational Warning: dedrooping was turned off');
        end
        
%Calculate electric field in sonde frame, according to 
%"physics convention" for positive E-component:
    %positive SE_x points from 1 to 0 (two lower channels)
    %positive SE_y points from 3 to 2 (two upper channels)
    %positive SE_z points upward
        
        SE_x=(V1_d-V0_d)*gamma/(2*G_xcal);
        SE_y=(V3_d-V2_d)*gamma/(2*G_ycal);
        SE_z=(V0_d+V1_d-V2_d-V3_d)*gamma/(4*G_zcal);
        
        SE_x_o=auto_offset(SE_x); %Get the data to start from zero
        SE_y_o=auto_offset(SE_y);
        SE_z_o=auto_offset(SE_z);
        
        Q_son=V0+V1+V2+V3;
                        
        figure(4);set(gcf,'color','w');


        if minute==20 
           if color==1
             a= plot(time_son,V0-0.03,'b-',time_son,V1-0.05,'k-',time_son,V2+0.03,'r-',time_son,V3+0.1,'m-',time_son,PPS/10,'k:'); 
             end;
           if color==0
             a= plot(time_son,V0-0.03,'k-',time_son,V1-0.05,'k-',time_son,V2+0.03,'k-',time_son,V3+0.1,'k-',time_son,PPS/10,'k:'); 
           end;
        end;
        xlim([ssec esec]);
        if minute==20 ylim([-0.6 0.53]); end
        title('Raw electrode voltages');
        ylabel('A/D Volts (V)');
        xlabel(['Time: Seconds since ',ZULUTIME,' UT']);
        set(a,'LineWidth',[1.5])
        if minute==20 set(gca,'XTick',[20.9 21.0 21.1  21.2 21.3 21.4 21.5 21.6 21.7 21.8 21.9 22.0]); end;
        if minute==20 str=['V_0' ;'V_1'; 'V_2';'V_3';'PPS'];x=[21.9 21.9 21.9 21.9 21.12]; y=[-0.09 -0.36 0.26 0.39 0.42]; end;
        h=text(x,y,str);
        loc='Southwest';

        Set_fonts(fontn,14,18,22,0,0,h,18);
        
        figure(5);set(gcf,'color','w');
        a=plot(time_son,SE_x_o,xc,time_son,SE_y_o+0.6,yc,time_son,SE_z_o-1,zc);
        set(a,'LineWidth',[1.5])
        xlim([ssec esec]);
        if minute==20 
            ylim([-4.5 1.5]); 
            xx=gca;
            set(xx,'XTick',[20.9 21.0 21.1 21.2 21.3 21.4 21.5 21.6 21.7 21.8 21.9 22.0])
        end
        title('Sonde-referenced E-field');
        ylabel('E-field (kV/m)');
        xlabel(['Time: Seconds since ',ZULUTIME,' UT']);
        if minute==20  x=[20.93 20.93 20.93]; y=[0.35 1.1 -0.35]; end;
       
        text(x(1),y(1),'E_{Sx}','FontName',fontn,'FontSize',18,'Color',xtc);
        text(x(2),y(2),'E_{Sy}','FontName',fontn,'FontSize',18,'Color',ytc);
        text(x(3),y(3),'E_{Sz}','FontName',fontn,'FontSize',18,'Color',ztc);
        
        h=0
        Set_fonts(fontn,14,18,22,0,0,h,18);
          
    %Calculate the electric field in Earth frame
        phi_T=phi_B+declination;     %Calculate angle from True North
        %phi_T=phi_true
        %All angles are still in degrees

        E_x=SE_x_o.*cosd(90-phi_T)-SE_y_o.*sind(90-phi_T); 
        E_y=SE_x_o.*sind(90-phi_T)+SE_y_o.*cosd(90-phi_T);
        E_z=SE_z_o; 
        %By using (90-phi_T) we are really doing a rotation matrix about an
        %angle that is positive counter-clockwise and starts at 0 on the
        %x-axis.
        
        figure(6);set(gcf,'color','w');  
        ymin=-0.5; ymax=2.5;

        if minute==20 xoff=3.3;yoff=0.3;zoff=2.5; end;
        if minute==12 xoff=1;yoff=0.2;zoff=0.5;end;
        a=plot(time_son,E_x+xoff,xc,time_son,E_y+yoff,yc,time_son,E_z+zoff,zc);
        set(a,'LineWidth',[1.5])
    %    xx=gca;
    %     set(xx,'XTick',[20.9 21.0 21.1 21.2 21.3 21.4 21.5 21.6 21.7 21.8 21.9 22.0])
        xlim([ssec esec]); 
         ylim([-1 4]); 
        
        title('Earth-referenced E-field');
        x=[20.93 20.93 20.93];
        y=[3.6 0.7 3.1];
        text(x(1),y(1),'E_X','FontName',fontn,'FontSize',18,'Color',xtc);
        text(x(2),y(2),'E_Y','FontName',fontn,'FontSize',18,'Color',ytc);
        text(x(3),y(3),'E_Z','FontName',fontn,'FontSize',18,'Color',ztc);
        h=0;
        %%%%
        str=['A' ;'B'; 'C';'D';'E';'F';'G';'H';'I'];
        x=[21.023 21.043 21.102 21.130 21.316 21.34 21.549 21.616 21.688]-0.013;
        x2=x+0.013;
        y=[2.02 1.92 1.68 1.5 1.31 1.25 0.1 -0.6 -0.84];
        h=text(x,y,str);
        y=[2.02 1.92 1.68 1.5 1.31 1.25 0.0 -0.7 -0.94];
        y2=y+0.25;
        y3=y2+0.25;
        
        line([x2 ;x2],[y2 ;y3],'LineWidth',[0.75],'Color',[0 0 0])
        %%%%
        
        ylabel('E-field (kV/m)');
        xlabel(['Time: Seconds since ',ZULUTIME,' UT']);
        Set_fonts(fontn,14,18,22,0,0,h,18);


        if draft==1 stampit; end

end

%************* atan2d ********************* 
function z = atan2d(y,x)
% atan2d.m --RS-function. Same as built-in function atan2, 
% but returns results in degrees instead of
% radians.
% Example
% atan2(1,1);
% 
% 45
z=(180.0/pi)*atan2(y,x);
end

%*************auto_offset ********************* 
function [output] = auto_offset( data)
%auto_offset -- 
%Assume that the first 1% of the input waveform is representative of the DC
%zero level.  Offsets rest of waveform by that amount.  Put in median for noise rejection.
sz=size(data);

if min(sz(1),sz(2))>1 
    error('You used an array that was not a simple vector.');
elseif max(sz(1),sz(2))==1
    error('You did not use an array.  You seem to have used a scalar');
end

    len=length(data);
    endindex=max(1,fix(0.01*len));  %Guarantee that last index is at least 1
    offset=median(data(1:endindex));
    output=data-offset;
end
   
%*************Set_fonts ********************* 
function a=Set_fonts(FigureFont,AxisFontSize,AxisLabelFontSize,...
    TitleFontSize,LegendHandle,LegendFontSize,TextHandle,TextFontSize)
%function a=Set_fonts(FigureFont,AxisFontSize,AxisLabelFontSize,...
%    TitleFontSize,LegendHandle,LegendFontSize,TextHandle,TextFontSize)
%
% FigureFont contains FontNames -- e.g. 'Times', 'Helevetica'
% Set's up the fonts of the last figure you plotted.  You should be done
% with the plotting and labeling before calling this function.
% Enter 0 for Handles you don't want to use. (TextHandle, LegendHandle)
L=LegendHandle;
T=TextHandle;
a=gca;
set(a,'FontName',FigureFont,'FontSize',AxisFontSize)
b=get(a,'Xlabel'); set(b,'FontName',FigureFont,'FontSize',AxisLabelFontSize)
b=get(a,'Ylabel'); set(b,'FontName',FigureFont,'FontSize',AxisLabelFontSize)
b=get(a,'Title'); set(b,'FontName',FigureFont,'FontSize',TitleFontSize)
if L~=0 set(L,'FontName',FigureFont,'FontSize',LegendFontSize); end
if T~=0  set(T,'FontName',FigureFont,'FontSize',TextFontSize); end
end

%*************lowpass ********************* 
function wave = lowpass( data,rate,tc )
%  s = lowpass( data,rate,tc ) e.g. s=lowpass(data,10000,1);
% Written by R.Sonnenfeld 11/2005
% Returns the FIR equivalent of an RC lowpass filter;
% Here data is an array of sampled data points (real laboratory data)
% rate is sample rate (10,000 is 10,000 samples/sec)
% tc is time constant of lowpass filter.  It is EXACTLY the value RC if it
% were a hardware RC filter.  Thus tc=0.1 means time const of 0.1 sec
% A lowpass filter is very simple:
%      a0=1-alpha==beta
%      b1= alpha
%      y_k=beta*x_k+alpha*y_k-1
% Here y_k is k-th sample of filtered data, x_k is k-th of input
% (unfiltered) data
% A highpass filter is only slightly more complex:
%      a0=(1+alpha)/2==gamma
%      a1=-1*gamma
%      b1= alpha
%      y_k=gamma*(x_k-x_k-1)+alpha*y_k-1
%
%      To find alpha: alpha=exp(-1/d) where d=rate*tc
%      For rate=10000 sample/sec, tc=0.1 sec, d=#samples to decay=1000
%      Since -1/1000<<1, then alpha=1-0.001=0.999, while 1-alpha=0.001
% Sample usage:
%   t=(0:10000)/10000;  Sample rate 10 kHz
%   y=sinewave(t,10,0); 10 Hz sinewave
%   plot(t,y,t,lowpass(y,10000,0.01)) Pass through filter with T=0.01 sec
d=rate*tc;
d_recip=-1/d; alpha=exp(d_recip); beta=1-alpha;
wave=data; %Create a wave as long as the input data array

% wave(1)=data(1); %Can't filter first point!
npts=length(data); 
for i=2:npts
    wave(i)=beta*data(i)+alpha*wave(i-1);
end
return 
end

%*************stampit ********************* 
function [fh,dh] = stampit(fd_flag)

% The Stampit function places a file and/or date stamp on the bottom
% of all open figures windows.  It must be called at the end of an m-file
% (otherwise the m-file name is []). If f_flag is 'f' then only the filename
% is placed on the figures, if d_flag = 'd' then only the date is stamped in
% the figure windows.  
%
% fh and dh are the corresponding handles of the filename and date text
% within each of the figure windows.  
%
% Call: [fh,dh]=stampit('f');
%
% Written by T. Burke, 8/21/01
%
% Revision List:
%  - addition of hidn_axis was added to properly locate the 
%    filename and date at the bottom of the overall figure window (DMP/TAB)
%  - addition of cur_ax in order to set the current axes handle back to its
%    original state (before Stampit is called) (DEK/TAB)
%  - addition of cur_fig in order to reset the current figure back to its 
%    original state (before Stampit was called) (TAB 11/7/01)
%  - place the handle for the hidden stampit axis to the bottom of the 
%    figure's Children handle list -- this should solve problems with the zoom
%    function and other functions after stampit has been called (TAB 11/13/01)
%  - used st(length(st)).name instead of st(2,1).name so that the original m-file
%    that is run is the one printed on the plots (ie. NOT a function's name in
%    case the figures are generated within a function call) (TAB 12/17/01)


%check for proper input to function
if nargin>1, error('Too Many Arguments for Stampit Function'); end
if nargin==0, flag=2; end
if nargin==1
    if fd_flag=='f'
        flag=1;
    else
        flag=0;
    end
end
 
%obtain location of calling m-file
[st,i] = dbstack;

if max(size(st))==1,error('Stampit must be called at the end of an m-file');end

h = findobj('Type','figure');
cur_fig = get(0,'CurrentFigure');

warning off
for i = 1:length(h)
    set(0,'CurrentFigure',h(i));
    cur_ax = get(h(i),'CurrentAxes'); % gets current axis handle for given figure
    
    hidn_axis=axes('Position',[0 0 1 1], 'Visible','off','Units','normalized','Tag','Stampit');
    % Make a hidden axis the size of the entire figure window
    
    if flag>0
        fh(i,1) = text(0.13,0.01,st(length(st)).name,'Interpreter','none'); % st(length(st)).name identifies the "parent" m-file called (TAB 12/17/01)
        set(fh(i,1),'FontSize',[9]);
    end    
    if flag==0|flag==2
        dh(i,1) = text(.905,0.01,['As of: ',date]);
        set(dh(i,1),'FontSize',[9],'HorizontalAlignment','Right','Interpreter','none');
    end
    
    set(h(i),'CurrentAxes',cur_ax); % Sets current axis back to its previous state
    
    chld = get(h(i),'Children');                              % These lines place the hidden "Stampit" axis
    if length(chld)>1                                         % to the bottom of the figure's Children handle list
        set(h(i),'Children',[chld(2:length(chld)); chld(1)])  % This alieveates problems with zoom and other functions
    end                                                       % after stampit has been called. (TAB 11/13/01)   
                                                          
end

set(0,'CurrentFigure',cur_fig);

pause(.001)
warning on
end