////////////////////////////////////////////////////////////////////////////////
/// name: schmidtplot                                                          /
/// filename: schmidtplot_1_0.sce                                              /
/// version: 1.0                                                               /
/// status: quick working version                                              /
/// author: lukas wischounig                                                   /
/// license: CeCILL-C                                                          /
/// date: march 2011                                                           /
/// platform: scilab 5.3.1 (written & tested on)                               /
/// os: win7-64 (written & tested on)                                          /
/// output: variables in the scilab workspace                                  /
/// keywords: tectonic, schmidt net, earth science, geology                    /
/// requires: matlab like plotting library                                     /
////////////////////////////////////////////////////////////////////////////////
///
/// description:
///
/// the script draws great circles and  points of a .xls file.
/// based on formulae after LISLE J.% R.LEYSHON (2004): Stereographic Projection Techniques for
/// Geologists and Civil Engineers - 2nd Edition. Cambridge University Press. 112 p.
///
/// for the basic structure of the input file have a look at the demo data files.
///
/// CONVENTION: 
/// +) tectonic data must follow CLAR convention
///
/// +) all angles must be in DEGREES!
///
/// +) the format of the input file must be
/// DipDir Dip Azimuth Plunge Color
///
/// where _DipDir / _Dip corresponds to GREAT CIRCLES
/// and _Azimuth / _Plunge corresponds to  POINTS
///
/// combinations of the two are also valid, ther result in the same color (or black = default)
///
/// the first line of the .xls file ("header") is IGNORED!!
/// xls files must be of type 97/2000/xp, i used openoffice
///
/// keep in mind that this version currently does not check input data 
/// for FORMAL CORRECTNESS. if input data are wrong you migth (or not) get an error 
/// i.e. a great circle at 387 34 will NOT cause an error! this means 387-360 = 27 34
/// have a look at demo_error_data if you are interested.....
///
/// color table
/// 1	blue
/// 2	green
/// 3	red
/// 4	cyan
/// 5	magenta
/// 6	brown/yellow
/// 7	gray
/// 8	black
///
/// TIP: to export the resulting plot to dxf i advise the use of gsview +pstoedit
/// (edit --> convert to vector format) after exporting the plot to eps or pdf. but don't be
/// disappointed, the hatch of the  point symbols has changed to a line in the
/// resulting dxf file.
///
////////////////////////////////////////////////////////////////////////////////
///
/// history:
///
/// version 1: quick programmed working version                       march 2011
///
////////////////////////////////////////////////////////////////////////////////

function [_rad_value] = _deg2rad(_deg_value)
    // angles(degrees) --> angles(radians)
    _rad_value = (%pi*_deg_value)/180;
endfunction

// some basic settings
_step = 100;
_radius = 1;
_counter = 0;

// get data
PathFileName_in = uigetfile(["*.xls"]);
[_path,_filename_in,_extension] = fileparts(PathFileName_in);

xls_sheets = readxls(PathFileName_in);
xls_sheet = xls_sheets(1);
[xls_s_heigth,xls_s_width]= size(xls_sheet)

// schmidt net: make a circle around the thing
_step_ci = 200;
_circ_x = cos(0:%pi/_step_ci:2*%pi); 
_circ_y = sin(0:%pi/_step_ci:2*%pi);
plot2d(_circ_x,_circ_y,8); // draw a circle around the whole thing

//  schmidt net: draw vertical & horizontal lines
_x_hor = [-1 1];
_y_hor = [0 0];
_x_vert = [0 0];
_y_vert = [-1 1.05]; 
plot2d(_x_hor,_y_hor,8);
plot2d(_x_vert,_y_vert,8);

_ha = gca(); //get handle of axes
_ha.isoview = 'on';

for _count = 2:xls_s_heigth // ignore line 1

    if ~isempty (xls_sheet(_count,1)) then // great circles
        _DipDir= %pi/2 - _deg2rad(xls_sheet(_count,1)); 
        _Dip= _deg2rad(xls_sheet(_count,2));

        if isempty(xls_sheet(_count,5)) then
            _color = 8;
        else
            _color =xls_sheet(_count,5);
        end
        // calculate great circles
        _D=[0:%pi/_step:%pi];
        _dp = _radius*sqrt(1-sin(atan(tan(_Dip).*sin(_D))));

        // plot great circles and rotate about (90°)-DipDir 
        plot2d(_dp.*sin(_D)*cos(_DipDir)-_dp.*cos(_D)*sin(_DipDir),...
        _dp.*sin(_D)*sin(_DipDir)+_dp.*cos(_D)*cos(_DipDir),_color);
        _counter = _counter +1;
    end

    if ~isempty (xls_sheet(_count,3)) then //  points
        _Azimuth = %pi/2 - _deg2rad(xls_sheet(_count,3)); 
        _Plunge = _deg2rad(xls_sheet(_count,4));

        if isempty(xls_sheet(_count,5)) then
            _color = 8;
        else
            _color =xls_sheet(_count,5);
        end
        // calculate  point
        _D=[%pi/2];
        _dp = _radius*sqrt(1-sin(atan(tan(_Plunge)*sin(_D))));

        // plot point and rotate about (90°)-_Azimuth 
        plot2d(_dp*sin(_D)*cos(_Azimuth)-_dp*cos(_D)*sin(_Azimuth),...
        _dp*sin(_D)*sin(_Azimuth)+_dp*cos(_D)*cos(_Azimuth));
        _e=gce();
        point=_e.children(1);
        point.mark_mode="on";
        point.mark_style=5;
        point.mark_foreground=8;
        point.mark_background=_color;
        _counter = _counter +1;
    end
end

// some plotting stuff, labels etc.
_ha = gca(); //get handle of axes
_ha.isoview = 'on';
_ha.axes_visible = ["off","off","off"]; // switch off ticks
_ti=_ha.title;
_ti.font_size = 3;
_ti.text="N";

_datastring = 'Datasets: '+ string(_counter); // datasets & filename labels
_datastring_2 = _filename_in;
xstring(0.35,-1,_datastring);
xstring(0.25,1,_datastring_2);










