DHOB (IU5SGN): Proiezioni stereografiche per geologia strutturale

mercoledì 11 luglio 2012

Proiezioni stereografiche per geologia strutturale

In questo post viene proposto un programma per disegnare le proiezioni stereografiche di utilizzo in geologia stratigrafica

Per prima cosa si deve calcolare l'intersezione del piano di cui si e' effettuata la misura (strato o discontinuita') con una sfera a centro (0,0,1) e raggio1; per fare cio' viene calcolato l'equazione del cerchio dell'equatore che poi viene ruotato dei corrispondenti valori di Strike e Dip
Dopo di cio' viene applicata la formula della proiezione stereografica per riportare i risultati sul piano di proiezione

----------------------
import pygame
import math

pygame.init()
screen = pygame.display.set_mode((300,300))

gamma = -math.pi # variabile per generare il cerchio massimo

teta = math.radians(90)   # angolo di dip
alfa = math.radians(159) # angolo di strike

raggio = 100
raggio2 = 50

pygame.draw.circle(screen,(255,0,0),(150,150),raggio,1)

while True:
    while gamma < (math.pi):
        #calcola il cerchio come con il centro in 0,0,1 e parallelo ad XY
        #in pratica calcola l'equatore come cerchio massimo
        x = math.cos(gamma)
        y = math.sin(gamma)
        z = 0
        #ruota il cerchio dell'equatore sull'asse x dell'angolo di pendenza (dip)
        x1 = x
        y1 = ((math.cos(teta)*y)-(math.sin(teta)*z))
        z1 = ((math.sin(teta)*y)+(math.cos(teta)*z))
        #ruota il cerchio dell'equatore del valore di azimuth
        x2 = ((math.cos(alfa)*x1)-(math.sin(alfa)*y1))
        y2 = ((math.sin(alfa)*x1)+(math.cos(alfa)*y1))
        z2 = 1+z1 #sposta l'origine della sfera a 0,0,1 per preparare la proiezione stereografica
        #proiezione stereografica
        if (z2<1): # per convenzione si plottano solo le intersezioni con la semisfera inferiore
            x3 = (2*x2)/(2-z2)
            y3 = (2*y2)/(2-z2)
            screen.set_at((int(x3*raggio2)+150,int(y3*raggio2)+150),(0,255,0))
        gamma = gamma+0.01
    pygame.display.flip()
----------------------