[GRASSLIST:4183] Re: structural geology - OFF TOPIC

Phillip J. Allen paallen at attglobal.net
Sun Aug 15 09:24:02 EDT 2004


Carlos,
I have written a function in pl/pgSQL (PostgreSQL) and Visual Basic (for 
Access) that takes Dip-Azimuth (0- <360 horizontal) and Vertical Dip 
Angle (0-90 downwards) and recalculates coordinates for a stereonet. 
 The sterionet coordinates have the origin 0,0 and a circle radius of 
100 unites.  So with your structural/orientation data in a database use 
this function in a query and open  the query in any graphing program and 
you have a stereonet.  I am not really sure if it is a wolff or schmitt 
net but I could never really tell  them appart anyway.  I always kept 
meaning to work these functions up into a R function.  

I have drawn the sterionet lines in Autocad and just open the dwg up in 
ArcView and can then sybolize the poles with any attribute I need.

Below are the PostgreSQL and Access functions.  Good luck.

I have always been meaning to work up a little script in 
PostgreSQL/PostGIS to create a 3D-Circle or Thin Disk representing the 
structural point and it's orientation ot view in 3-d but I have just 
never gotten around to it yet.  

Good Luck.


Phillip J. Allen
Chief Geochemist Anglogoldashanti South America
Lima Peru
e-mail: paallen at attglobal.net


Carlos Henrique Grohmann wrote:

>Hello, sorry about the off-topic, but dows anyone kkow a software for structural
>geology (stereonet) for linux? I have beeen looking on the net, but haven't
>find nothing yet.
>
>thanks
>
>Carlos
>  
>
****************************
*****************************
**PostgreSQL PL/PGSQL Function**
*****************************
*****************************
-- Name: "stereonet" (smallint,double precision,double precision) Type: 
FUNCTION Owner: paallen
--

CREATE FUNCTION "stereonet" (smallint,double precision,double precision) 
RETURNS double precision AS '
       DECLARE
           xyz ALIAS FOR $1;
           dipazimuth FLOAT8;
           dipangle FLOAT8;
           originx CONSTANT FLOAT8 := 0.0;
           originy CONSTANT FLOAT8 := 0.0;
           circdiam CONSTANT FLOAT8 := 100.0;
           spherex FLOAT8;
           spherey FLOAT8;
           spherez FLOAT8;
           stnetx FLOAT8;
           stnety FLOAT8;
       BEGIN
-- Calculate XY Coordinates for Wulff Steronet with origin at 0,0,0 and 
circle diameter of 100
--           Convert Plane to Pole
           dipazimuth = $2;
           dipangle = $3;
           dipazimuth = dipazimuth + 180.0;
           IF dipazimuth >= 360.0 THEN
          dipazimuth = dipazimuth  -360.0;
           END IF;
           dipangle = 0 - (90.0 - dipangle);
          
           spherex = returnpc2newxyz(originx, 
originy,0.0,circdiam,dipazimuth, dipangle,1);
           spherey = returnpc2newxyz(originx, 
originy,0.0,circdiam,dipazimuth,dipangle,2);
           spherez = returnpc2newxyz(originx, 
originy,0.0,circdiam,dipazimuth,dipangle,3);
          
           IF spherez >=0 THEN
          stnetx = spherex;
          stnety = spherey;
           ELSE
              stnetx = (((0 - circdiam) / (spherez - 100)) * (spherex - 
originx)) + originx;
          stnety = (((0 - circdiam) / (spherez - 100)) * (spherey - 
originy)) + originy;
           END IF;

           IF xyz = 1 THEN
          RETURN stnetx;
           ELSIF xyz = 2 THEN
          RETURN stnety;
           ELSIF xyz = 3 THEN
              RETURN 0.0;
           ELSIF xyz = 4 THEN
          RETURN spherex;
           ELSIF xyz = 5 THEN
          RETURN spherey;
           ELSIF xyz = 6 THEN
              RETURN spherez;
           ELSE
              RETURN -9999999.0;
           END IF;
     END;
' LANGUAGE 'plpgsql';

--
-- Name: "returnpc2newxyz" (double precision,double precision,double 
precision,double precision,double precision,double precision,smallint) 
Type: FUNCTION Owner: paallen
--

CREATE FUNCTION "returnpc2newxyz" (double precision,double 
precision,double precision,double precision,double precision,double 
precision,smallint) RETURNS double precision AS '
      DECLARE
        newx FLOAT8;
        newy FLOAT8;
        newz FLOAT8;
        xcenter ALIAS FOR $1;
        ycenter ALIAS FOR $2;
        zcenter ALIAS FOR $3;
        dist ALIAS FOR $4;
        azimuthdd ALIAS FOR $5;
        inclinangdd ALIAS FOR $6;
        horizdist FLOAT8;
        ddrad CONSTANT FLOAT8 := 3.14159265358979/180.0;
       BEGIN
      
        IF inclinangdd < 0 THEN
           newz = zcenter - (dist * sin((abs(inclinangdd)) * ddrad));
        ELSE
           newz = zcenter + (dist * sin((abs(inclinangdd)) * ddrad));
        END IF;

        IF $7 = 3 THEN
               RETURN newz;
        END IF;
       
        horizdist = dist * cos(abs(inclinangdd * ddrad));
       
        IF $7 = 1 THEN
           newx = xcenter + (horizdist * (sin(azimuthdd * ddrad)));
           RETURN newx;
        ELSIF $7 = 2 THEN
           newy = ycenter + (horizdist * (cos(azimuthdd * ddrad)));
           RETURN newy;
        END IF;
        RETURN -9999999.0;
    END;
' LANGUAGE 'plpgsql';






*****************************
*****************************
**Access Visual Basic Functions**
*****************************
*****************************

Function StereoNet(ByVal snt As String, xyz As String, dazimuth As 
Double, dangle As Double) As Double
Rem*************************************************************************************************************
Rem*********Calculate Wulff or Lambert X Stereo Net Coordinates from Dip 
Azimuth & Dip Angle degrees************
Rem*************************************************************************************************************
Rem****************************Great Circle is 100 units 
Radius*************************************************

Dim originx, originy, circdiam As Double
Dim spherex, spherey, spherez As Double
Dim stnetx, stnety As Double
Dim anb As Integer

originx = 0
originy = 0
circdiam = 100

Rem Convert plane to pole
dazimuth = dazimuth + 180
If (dazimuth >= 360) Then
    dazimuth = dazimuth - 360
End If

dangle = 90 - dangle

Rem anb = MsgBox(dazimuth, vbOKOnly, "azimuth")
Rem anb = MsgBox(dangle, vbOKOnly, "angle")


Rem If (snt = "W") Then Rem Wulff Stereo Net
    spherex = ReturnPC2NewX(originx, originy, 0, circdiam, dazimuth, 0 - 
dangle)
    spherey = ReturnPC2NewY(originx, originy, 0, circdiam, dazimuth, 0 - 
dangle)
    spherez = ReturnPC2NewZ(originx, originy, 0, circdiam, dazimuth, 0 - 
dangle)
Rem anb = MsgBox(spherex, vbOKOnly, "spherex")
Rem anb = MsgBox(spherey, vbOKOnly, "spherey")
Rem anb = MsgBox(spherez, vbOKOnly, "spherez")
   
    If (spherez >= 0) Then
   
        stnetx = spherex
        stnety = spherey
    Else
        stnetx = (((0 - circdiam) / (spherez - 100)) * (spherex - 
originx)) + originx
        stnety = (((0 - circdiam) / (spherez - 100)) * (spherey - 
originy)) + originy
       
    End If
   
Rem ElseIf (snt = "L") Then Rem Lambert Stereo Netdssd


Rem End If

If xyz = "x" Then
    StereoNet = stnetx
ElseIf xyz = "y" Then
    StereoNet = stnety
ElseIf xyz = "sx" Then
    StereoNet = spherex
ElseIf xyz = "sy" Then
    StereoNet = spherey
ElseIf xyz = "sz" Then
    StereoNet = spherez
End If


End Function


Public Function ReturnPC2NewY(ByVal XCPoint, YCPoint, ZCPoint, dist, 
AzimuthAnglDD, InclinAnglDD)
Dim NewZ As Double
Dim DH As Double

Rem Calculate Depth(z)+z is upward, 0 is horizontal, -z is downward
If InclinAnglDD < 0 Then
    NewZ = ZCPoint - (dist * Sin((Abs(InclinAnglDD)) * 3.14159265358979 
/ 180))
Else
    NewZ = ZCPoint + (dist * Sin((Abs(InclinAnglDD)) * 3.14159265358979 
/ 180))
End If

Rem Calculate Horizontal X Distance
DH = dist * Cos(Abs(InclinAnglDD * 3.14159265358979 / 180))

ReturnPC2NewY = YCPoint + (DH * (Cos(AzimuthAnglDD * 3.14159265358979 / 
180)))
End Function

Function ReturnPC2NewZ(ByVal XCPoint, YCPoint, ZCPoint, dist, 
AzimuthAnglDD, InclinAnglDD)
Dim DH As Double

Rem Calculate Depth(z)+z is upward, 0 is horizontal, -z is downward
If InclinAnglDD < 0 Then
    ReturnPC2NewZ = ZCPoint - (dist * Sin((Abs(InclinAnglDD)) * 
3.14159265358979 / 180))
Else
    ReturnPC2NewZ = ZCPoint + (dist * Sin((Abs(InclinAnglDD)) * 
3.14159265358979 / 180))
End If

End Function

Function ReturnPC2NewX(ByVal XCPoint, YCPoint, ZCPoint, dist, 
AzimuthAnglDD, InclinAnglDD)
Dim NewX As Double, NewZ As Double
Dim DH As Double

Rem Calculate Depth(z)+z is upward, 0 is horizontal, -z is downward
If InclinAnglDD < 0 Then
    NewZ = ZCPoint - (dist * Sin((Abs(InclinAnglDD)) * 3.14159265358979 
/ 180))
Else
    NewZ = ZCPoint + (dist * Sin((Abs(InclinAnglDD)) * 3.14159265358979 
/ 180))
End If

Rem Calculate Horizontal X Distance
DH = dist * Cos(Abs(InclinAnglDD * 3.14159265358979 / 180))

ReturnPC2NewX = XCPoint + (DH * (Sin(AzimuthAnglDD * 3.14159265358979 / 
180)))
End Function






More information about the grass-user mailing list