stereo_plane.m

function stereo_plane(Attitude,Dip,DipD)
%Plots the stereographic projection of a plane using a Wulff stereonet.
%   stereo_plane(Attitude,Dip,DipD) plots a plane over a Wulf stereonet
%   and evaluates the direction of the dip angle in order to choose the
%   right major circle. Attitude and Dip are both on sexagecimal degrees
%   while DipD is a logical parameter that can only be 'N', 'S', 'E' or 'W'
%   and single quotes shall be respected. Attitude, Dip and DipD can be
%   vectors in order to evaluate different planes on a single figure, but
%   these vectors have to be the same lenght, in case the plane is at 0,180
%   or 360, DipD must be 'E' or 'W'.
%
%   If dip direction is known, you can use the function included 
%   on this package called dipd2str in order to get the alphabetic code
%   for DipD

    North=cell(length(DipD),1);
    East=cell(length(DipD),1);
    legh=zeros(length(DipD),1);
    message2=cell(length(DipD),1);
%    legvar=cell(length(DipD),1);
    for i=1:length(DipD)
        if DipD(i)=='N'
           North{i}=true;
        elseif DipD(i)=='S'
           North{i}=false;
        elseif DipD(i)=='E'
           East{i}=true;
        elseif DipD(i)=='W'
           East{i}=false;
        end
    end
    figure
    if length(Attitude)==1
        ftitle='Plane over the Wulff projection';
    else
        ftitle='Planes over the Wulff projection';
    end
    set(gcf,'units','normalized','position',[0 0 1 1],'name',ftitle,'NumberTitle','off')
    No = 50;                                           
    xh = [-sind(Attitude) sind(Attitude)];yh = [-cosd(Attitude) cosd(Attitude)];
    axis([-1 1 -1 1]);axis('square');
    rectangle('Position',[-1 -1 2 2],'Curvature',[1 1],'FaceColor','w')            
    axis off;hold on;
    if Dip==90
        plot(xh,yh,'-k');
    else
        psi = 0:180/No:180;%Puntos del circulo mayor, de cero a 180 grados
        for i=1:length(Dip)
            rdip = Dip(i);                          
            radip = atand(tand(rdip).*sind(psi));
            rproj = tand((180/2 - radip)/2);
%            rproj0 = tand((180/2 - atand(tand(rdip).*sind(0)))/2);
            x1 = rproj .* sind(psi+Attitude(i));
            x2 = rproj .* (-sind(psi-Attitude(i)));
            y1 = rproj .* cosd(psi+Attitude(i));
            y2 = rproj .* cosd(psi-Attitude(i));
            if Attitude(i)==0||Attitude(i)==360
                if East{i}==true
                    legh(i)=plot(x1,y1,'-k');
                elseif East{i}==false
                    legh(i)=plot(x2,y2,'-k');
                end
            elseif Attitude(i)==180
                if East{i}==true
                    legh(i)=plot(x2,y2,'-k');
                elseif East{i}==false
                    legh(i)=plot(x1,y1,'-k');
                end
            else
                if rproj * cosd(180/2+Attitude(i))>rproj * cosd(180/2-Attitude(i));
                    %p1 is N p2 is S
                    if North{i}==true
                        legh(i)=plot(x1,y1,'-k');
                    elseif North{i}==false
                        legh(i)=plot(x2,y2,'-k');
                    end
                elseif rproj * cosd(180/2+Attitude(i))<rproj * cosd(180/2-Attitude(i))
                        %p2 is N p1 is S
                    if North{i}==true
                        legh(i)=plot(x2,y2,'-k');
                    elseif North{i}==false
                        legh(i)=plot(x1,y1,'-k');
                    end
                end
            end            
            message=sprintf('(%d) Strike: %.2fº Dip: %.2fº%c\n',i,Attitude(i),Dip(i),DipD(i));
            message2{i}=horzcat(message);
        end
%       legvar
        %leg=legend(legh,message2);
        %set(leg,'location','northeastoutside')
    end
    coordinates={'N';'E';'S';'W'};verthor={'VerticalAlignment';'HorizontalAlignment';'VerticalAlignment';'HorizontalAlignment'};
    al={'Bottom';'Left';'Top';'Right'};xtext= [0 1 0 -1];ytext= [1 0 -1 0];
    for i=1:4
        text(xtext(i),ytext(i),coordinates(i),verthor(i),al(i))
    end
    text(-2.4,0,message2(1:end/2),'HorizontalAlignment','Left','BackgroundColor','w','FontSize',8)
    text(-1.7,0,message2(end/2+1:end),'HorizontalAlignment','Left','BackgroundColor','w','FontSize',8)
end