%**************************************************************************
% CLUSTER project: Plot divB & curlB to screen     Autor: Christian Maszl *
% Unit System = [SI]                  Mail: e0025754@student.tuwien.ac.at *
%**************************************************************************

% This subroutine plot jav and B on the screen. The program looks
% up the datafiles in the "outputdata" directrory defined in
% config.m . Please make sure that you selected the right mode!

% Load configfile and check the path's ...
syscheck;
date=input('Please enter the desired day [yymmdd]: ','s');
% Extract extension
l=size(datatype,2);
if (l==4)
    exten=datatype(2:4);
else
    exten=datatype(2:4);
end
javname =[outputspath date '0000' mode '.jav.gse.hr.igm'  exten];
B14name =[outputspath date '0000' mode '.B1-4.gse.hr.igm' exten];
divBname=[outputspath date '0000' mode '.divB.gse.hr.igm' exten];
% Load the variabels to the memory
try
load(javname,datatype);
    disp('jav loaded to memory!');
load(B14name,datatype);
    disp('B1-B4 loaded to memory!');
    B1=B(:,2:4);B2=B(:,5:7);B3=B(:,8:10);B4=B(:,11:13);
load(divBname,datatype);
    disp('<divB> loaded to memory!');
catch
    disp('Error: Filename not found! Please check your input or the available files!');
    disp('Please restart the program!');
    break;
end
clear javname B14name divBname l exten year month day;
% Convert the date for display use
try
year(1,1) =str2num(date(1,1:2));
month(1,1)=str2num(date(1,3:4));
day(1,1)  =str2num(date(1,5:6));
datum=datestr(datenum([2000+year month day]));
T=jav(:,1);
% Find plot range for jav
range =(sum(abs(jav) ,1)/size(jav,1))*10;
% Define styles, labels for axis, ...
jav_axisr=[0 24.0 -range(1,2) range(1,2)];
timelabel=['{t / [h]}'];
% Plot the x component of the current density versus the time.
subplot(4,1,1)
plot(T,jav(:,2))
    axis(jav_axisr)
    ylabel('{J}_{x}{ / [A/m}^{2}{]}','FontSize',10)
    xlabel(timelabel    ,'FontSize',10)
    text1=[datum ' , Coord. system: GSE , Unit System: SI , Time System:' ...
		' UT'];    
    title(text1,'FontSize',12)
    grid on
% Plot the y component of the current density versus the time.
subplot(4,1,2)
plot(T,jav(:,3)) 
    axis(jav_axisr)
    ylabel('{J}_{y}{ / [A/m}^{2}{]}','FontSize',10)
    xlabel(timelabel    ,'FontSize',10)
    grid on
% Plot the z component of the current density versus the time.
subplot(4,1,3)
plot(T,jav(:,4)) 
    axis(jav_axisr)
    ylabel('{J}_{z}{ / [A/m}^{2}{]}','FontSize',10)
    xlabel(timelabel    ,'FontSize',10)
    grid on
% Plot <divB>/abs(curlB) versus the time.
divBocurlB=(divB(:,2)./sqrt(dot(jav(:,2:4),jav(:,2:4),2)))*100;
range =(sum(abs(divBocurlB) ,1)/size(jav,1))*10;
subplot(4,1,4)
plot(T,divBocurlB) 
    axis([0.0 24.0 0 range]) 
    ylabel('{\nabla\cdotB/|\nablaxB| / [%]}' ,'FontSize',10)
    xlabel(timelabel            ,'FontSize',10)
    formula1=['{\it\mu}_{\it0}{\itJ}_{\itav}{=}{\it\nablaxB}'];
    formula2=['{\it\int_V}{\itdV \nabla\cdotB}={\it\int_A}{\itdA B}'];
    text=['Used formula:   ' formula1 '   and   ' formula2 '     Data' ...
         ' source: {\itFGM}'];
    title(text,'Fontsize',12)
    grid on
% Open a new plot window
figure
% Calculate average of the B fiels
BX=(B1+B2+B3+B4)/4;
% Plot the x component of B versus the time
subplot(4,1,1)
plot(T,BX(:,1)) 
    axis([0.0 24.0 -1000 1000])
    axis 'auto y'
    ylabel('{B}_{x}{ / [nT]}','FontSize',10)
    xlabel(timelabel  ,'FontSize',10)
    title(text1,'FontSize',12)
    grid on
% Plot the y component of B versus the time.
subplot(4,1,2)
plot(T,BX(:,2)) 
    axis([0.0 24.0 -1000 1000])
    axis 'auto y'
    ylabel('{B}_{y}{ / [nT]}','FontSize',10)
    xlabel(timelabel  ,'FontSize',10)
    grid on
% Plot the z component of B versus the time.
subplot(4,1,3)
plot(T,BX(:,3)) 
    axis([0.0 24.0 -1000 1000])
    axis 'auto y'
    ylabel('{B}_{z}{ / [nT]}','FontSize',10)
    xlabel(timelabel  ,'FontSize',10)
    grid on
% Calculate |B|
B=sqrt(dot(BX,BX,2));
% Plot |B| versus the time.
subplot(4,1,4)
plot(T,B) 
    axis([0.0 24.0 -1000 1000])
    axis 'auto y'
    ylabel('{|B| / [nT]}' ,'FontSize',10)
    xlabel(timelabel    ,'FontSize',10)
    text=['Used formula:   ' '{\itB=(B}_{\it1}{+\itB}_{\it2}+'...
	  '{\itB}_{\it3}{+\itB}_{\it4}{\it)/4}' '    Data source:  '...
	  '{\itFGM}'];
    title(text,'FontSize',12)
    grid on
catch
   disp('Error: Files contain no data!')
end
clear
