%**************************************************************************
% CLUSTER project: Cluster Flight                  Autor: Christian Maszl * 
% Unit System = [SI]                  Mail: e0025754@student.tuwien.ac.at *
%**************************************************************************
%
% This program simulates flights of the four cluster spacecrafts
% trough various current sheets.

% Define constants
y0 = 1.256637061e-6; %[y0]=Vs/Am
disp ('Please choose the type of the current sheet: ');
message=['[1] Cylindrical current sheet - [2] Flat current sheet'];
disp (message);
sheetform=input('Please type [1] or [2]: ');
disp(' ');
if (sheetform==1)
   r =input('Radius of the current tube [m]: r= ');
   s =input('Satellite separation       [m]: s= ');
   jy=input('Current density          [A/m²]: j= ');
   disp(' ');
   s=s/2;
% Compute the x,y and z values depending on the flight path.
        [x,y,z,T,stepanz]=bewegungsrichtung(r,s);
% Compute the position vectors of the space crafts.
        [r1,r2,r3,r4]=sat_pos(x,y,z,s);
% Report when the spacecrafts enter the current sheet
%                      enter_sheet(r1,r2,r3,r4,r);
% Compute the B-field vectors at the tetrahedron vertices.
        [B1]=cylsheet_B(r1,r,y0,jy,stepanz);
        [B2]=cylsheet_B(r2,r,y0,jy,stepanz);
        [B3]=cylsheet_B(r3,r,y0,jy,stepanz);
        [B4]=cylsheet_B(r4,r,y0,jy,stepanz);
% Use the curlometer method to compute jy again.
        r1=r1*1e-3;r2=r2*1e-3;r3=r3*1e-3;r4=r4*1e-3;
        [jav,divB]=curlB(B1,B2,B3,B4,r1,r2,r3,r4,y0);
% Plot the results on the screen
        jav=[T jav];
        plotdiagram(r1,r2,r3,r4,B1,B2,B3,B4,jav,divB,stepanz);
elseif (sheetform==2)
   d =input('Thickness of the flat current sheet [m]: d= ');
   s =input('Satellite separation                [m]: s= ');
   jy=input('Current density                   [A/m²]: j= ');
   s=s/2; d=d/2;
% Compute the x,y and z values depending on the flight path.
        [x,y,z,T,stepanz]=bewegungsrichtung(d,s);
% Compute the position vectors of the space crafts.
        [r1,r2,r3,r4]=sat_pos(x,y,z,s);
% Compute the B-field vectors at the tetrahedron vertices.
        [B1]=flatsheet_B(r1,d,y0,jy,stepanz);
        [B2]=flatsheet_B(r2,d,y0,jy,stepanz);
        [B3]=flatsheet_B(r3,d,y0,jy,stepanz);
        [B4]=flatsheet_B(r4,d,y0,jy,stepanz);
% Use the curlometer method to compute jy again.
        r1=r1*1e-3;r2=r2*1e-3;r3=r3*1e-3;r4=r4*1e-3;
        [jav,divB]=curlB(B1,B2,B3,B4,r1,r2,r3,r4,y0);
% Plot the results on the screen
        r1=r1;r2=r2;r3=r3;r4=r4;
        jav=[T jav];
        plotdiagram(r1,r2,r3,r4,B1,B2,B3,B4,jav,divB,stepanz);
        r1=r1*1e3;r2=r2*1e3;r3=r3*1e3;r4=r4*1e3;
else  
end
