Calcolo numerico : sistema terrestre locale

E-school di Arrigo Amadori

Calcolo numerico

Sistema terrestre locale

Per maggiori informazioni :

http://www.arrigoamadori.com/lezioni/Miscellanea/1/SistemaTerrestreLocale.htm
http://www.arrigoamadori.com/lezioni/Miscellanea/1/Guglielmini.htm .

parametri del sistema :

T = (periodo di rotazione della Terra)
G = (costante di gravitazione universale)
M = (massa della Terra)
R0 = (raggio della Terra)
θ = (colatitudine in radianti) (*)
φ = (longitudine in radianti)
ro = (densità dell'atmosfera solidale con il sistema ruotante)
k1 = (coeff. attrito del mezzo (per parte dipendente dalla velocità lineare))
k2 = (coeff. attrito del mezzo (per parte dipendente dalla velocità al quadrato))

parametri iniziali particella sferica rispetto al sistema terrestre locale :

	massa :	raggio :	posizione :			velocità :
			x	y	z	x	y	z
particella :

parametri di iterazione :

dt = (incremento tempo) ni = (numero iterazioni)

parametri di scala :

lar = alt = (larghezza ed altezza immagine in pixel)

xa = 0         xb = (estremi asse x)
ya = 0         yb = (estremi asse y)
za = 0         zb = (estremi asse z)

Operatori ammessi : - (per numeri negativi) . (per numeri decimali)

(*) i dati di defaults (eccetto il raggio del grave) sono quelli dell'esperimento di G.B. Guglielmini (1790, Specola bolognese)

Sorgenti PHP :

sistematerrlocale.php

<html>

<head>
<title>Sistema terrestre locale</title>
</head>

<body>

<script language="php">

error_reporting (E_ALL ^ E_NOTICE);

// parametri in input
//
// $T = periodo rotazione della Terra
// $G = costante di gravitazione universale
// $M = massa della Terra
// $R0 = raggio della Terra
// $teta0 = colatitudine in radianti
// $fi0 = longutudine in radianti
// $ro = densità atmosfera
// $k1 = coeff. attrito del mezzo (per parte dipendente dalla velocità lineare)
// $k2 = coeff. attrito del mezzo (per parte dipendente dalla velocità al quadrato)
// particella : $x0 = posizione iniziale x $y0 = posizione iniziale y $z0 = posizione iniziale z
// particella : $vx0 = velocità iniziale x $vy0 = velocità iniziale y $vz0 = velocità iniziale z
// $m = massa particella
// $r = raggio particella
// $lar = larghezza immgine in pixel
// $alt = altezza immagine in pixel
// $xb = fine dominio x
// $yb = fine dominio y
// $zb = fine dominio z
// $ni = numero iterazioni
// $dt = incremento tempo

// echo $G . " " . $R0 . " " . $M . "";

//
// parametri in output
//
// crea l'immagine PNG e ne visualizza il link
//

// inserisce tutti gli include (funzioni esterne)

if (is_integer(strpos($SERVER_SOFTWARE, "Win")))
{
//include("..\Routines\derivataParz.php");
//include("..\Routines\daPowAVb.php");
//include("..\Routines\daVbAPow.php");
//include("..\Routines\derivataMonom.php");
//include("..\Routines\semplEspr.php");
}
else
{
//include("../Routines/derivataParz.php");
//include("../Routines/daPowAVb.php");
//include("../Routines/daVbAPow.php");
//include("../Routines/derivataMonom.php");
//include("../Routines/semplEspr.php");
}

// controlli vari

// parametri di scala

$xa = 0; // inizio dominio x
$ya = 0; // inizio dominio y
$za = 0; // inizio dominio z

$dx = ($xb - $xa) / ($lar / 2 - 1); // intervallo di scomposizione asse x
$dy = ($yb - $ya) / ($lar / 2 - 1); // intervallo di scomposizione asse y
$dz = ($zb - $za) / ($alt / 2 - 1); // intervallo di scomposizione asse z

$alfa = atan($lar / $alt); // angolo fra asse x e parte negativa asse z
$coeff = sqrt($lar*$lar/4+$alt*$alt/4)/$xb; // coefficiente prospettico

// disegno frame con coordinate

$im = ImageCreate ($lar + 100, $alt + 100);

$bianco = ImageColorAllocate ($im, 255, 255, 255);
$nero = ImageColorAllocate ($im, 0, 0, 0);
$rosso = ImageColorAllocate ($im, 255, 0, 0);
$verde = ImageColorAllocate ($im, 0, 255, 0);
$blu = ImageColorAllocate ($im, 0, 0, 255);

Imageline ($im, 0, 0, $lar, 0, $bianco);
Imageline ($im, 0, $alt, $lar, $alt, $bianco);
Imageline ($im, 0, 0, 0, $alt, $bianco);
Imageline ($im, $lar, 0, $lar, $alt, $bianco);

Imageline ($im, $lar / 2, 0, $lar / 2, $alt / 2, $nero);
Imageline ($im, $lar / 2, $alt / 2, $lar, $alt / 2, $nero);
Imageline ($im, $lar / 2, $alt / 2, 0, $alt, $nero);

ImageString ($im, 1, $lar / 2, $alt / 2 + 10, $xa, $nero);
ImageString ($im, 1, 0, $alt + 10, $xb . " (x)" , $nero);
ImageString ($im, 1, $lar / 2 + 10, $lar / 2 - 10, $ya, $nero);
ImageString ($im, 1, $lar + 10, $alt / 2 - 10, $yb . " (y)", $nero);
ImageString ($im, 1, $lar / 2 - 10, $lar / 2 - 10, $za, $nero);
ImageString ($im, 1, $lar / 2 - 10, 0, $zb . " (z)", $nero);

ImageString ($im, 2, 0, $alt + 30, "T = " . $T . " G = " . $G . " M = " . $M . " R0 = " . $R0 . " teta = " . $teta0 . " fi = " . $fi0, $nero);

ImageString ($im, 2, 0, $alt + 50, "ro = " . $ro . " k1 = " . $k1 . " k2 = " . $k2 . " m = " . $m . " r = " . $r, $nero);

ImageString ($im, 2, 0, $alt + 70, "Valori iniziali : x = " . $x0 . " y = " . $y0 . " z = " . $z0 . " vx = " . $vx0 . " vy = " . $vy0 . " vz = " . $vz0 . " R = " . sqrt($x0*$x0+$y0*$y0+$z0*$z0) . " v = " . sqrt($vx0*$vx0+$vy0*$vy0+$vz0*$vz0), $nero);

$omega = 2*M_PI/$T;

$cost = $G*$M*((4/3)*M_PI*$r*$r*$r*$ro/$m - 1);

$g = $G*$M/($R0*$R0); // acelerazione di gravità locale
$deltax = $z0*sin($teta0)*cos($teta0)*pow($omega,2)*$R0/sqrt(pow($omega,4)*pow($R0,2)*pow(sin($teta0),2)-2*$g*pow($omega,2)*$R0*pow(sin($teta0),2)+pow($g,2)); // deviazione del filo a piombo sulla verticale

echo "--- dati di riferimento validi solo per la semplice caduta di un grave :";

echo "g = " . $g . " (accelerazione di gravità al suolo)";

echo " deviazione teorica verso sud del filo a piombo = " . $deltax . " (inizio filo in (0,0,z) lunghezza filo = z)";

$devteor = (1/3)*pow(2*$z0/$g,3/2)*$g*$omega*cos(M_PI/2-$teta0);

echo "deviazione teorica verso est secondo Landau = " . $devteor . " (caduta di un grave da quota z senza atmosfera e trascurando la forza centrifuga)";

//echo $cost . "";

$fx = "$m*($cost*(\$x)/pow((\$x)*(\$x)+(\$y)*(\$y)+(\$z)*(\$z),3/2)-((\$vx)/$m)*($k1+$k2*sqrt((\$vx)*(\$vx)+(\$vy)*(\$vy)+(\$vz)*(\$vz)))+2*$omega*(\$vy)+$omega*$omega*(\$x))";

$fy = "$m*($cost*(\$y)/pow((\$x)*(\$x)+(\$y)*(\$y)+(\$z)*(\$z),3/2)-((\$vy)/$m)*($k1+$k2*sqrt((\$vx)*(\$vx)+(\$vy)*(\$vy)+(\$vz)*(\$vz)))-2*$omega*(\$vx)+$omega*$omega*(\$y))";

$fz = "$m*($cost*(\$z)/pow((\$x)*(\$x)+(\$y)*(\$y)+(\$z)*(\$z),3/2)-((\$vz)/$m)*($k1+$k2*sqrt((\$vx)*(\$vx)+(\$vy)*(\$vy)+(\$vz)*(\$vz))))";

//echo $fx . "";
//echo $fy . "";
//echo $fz . "";

// settaggio matrici valori iniziali

$ex1 = cos($teta0) * cos($fi0);
$ex2 = cos($teta0) * sin($fi0);
$ex3 = -sin($teta0);

$ey1 = (1 / abs(sin($teta0))) * (-1) * sin($teta0) * sin($fi0);
$ey2 = (1 / abs(sin($teta0))) * sin($teta0) * cos($fi0);
$ey3 = 0;

$ez1 = sin($teta0) * cos($fi0);
$ez2 = sin($teta0) * sin($fi0);
$ez3 = cos($teta0);

$R0x = $R0 * sin($teta0) * cos($fi0);
$R0y = $R0 * sin($teta0) * sin($fi0);
$R0z = $R0 * cos($teta0);

$Rx = $R0x + $x0 * $ex1 + $y0 * $ey1 + $z0 * $ez1;
$Ry = $R0y + $x0 * $ex2 + $y0 * $ey2 + $z0 * $ez2;
$Rz = $R0z + $x0 * $ex3 + $y0 * $ey3 + $z0 * $ez3;

$Vx = $vx0 * $ex1 + $vy0 * $ey1 + $vz0 * $ez1;
$Vy = $vx0 * $ex2 + $vy0 * $ey2 + $vz0 * $ez2;
$Vz = $vx0 * $ex3 + $vy0 * $ey3 + $vz0 * $ez3;

//echo $R0x . " " . $R0y . " " . $R0z . "";

//echo $Rx . " " . $Ry . " " . $Rz . "";

//$u0["x"] = $x0; $u0["y"] = $y0; $u0["z"] = $z0; // posizioni
$u0["x"] = $Rx; $u0["y"] = $Ry; $u0["z"] = $Rz; // posizioni

//$v0["x"] = $vx0; $v0["y"] = $vy0; $v0["z"] = $vz0; // velocità
$v0["x"] = $Vx; $v0["y"] = $Vy; $v0["z"] = $Vz; // velocità

$t = 0;
for ($ii = 1; $ii <= $ni; $ii++) // iterazione
{
$t = $t + $dt;

// disegna posizione particella

$x = $u0["x"];
$y = $u0["y"];
$z = $u0["z"];

$xprimo = $x - $R0x;
$yprimo = $y - $R0y;
$zprimo = $z - $R0z;

$xxprimo = $ex1 * $xprimo + $ex2 * $yprimo + $ex3 * $zprimo;
$yyprimo = $ey1 * $xprimo + $ey2 * $yprimo + $ey3 * $zprimo;
$zzprimo = $ez1 * $xprimo + $ez2 * $yprimo + $ez3 * $zprimo;

//echo $xxprimo . " " . $yyprimo . " " . $zzprimo . "";

if ($xxprimo >= -$xb and $xxprimo <= $xb and $yyprimo >= -$yb and $yyprimo <= $yb and $zzprimo >= -$zb and $zzprimo <= $zb)
{
$xx = $xxprimo - $xa;
$yy = $lar / 2 + (($yyprimo - $ya)) / $dy;
$zz = $alt / 2 - (($zzprimo - $za)) / $dz;
$yyy = $yy - $coeff * $xx * sin($alfa);
$zzz = $zz + $coeff * $xx * cos($alfa);
Imageline ($im, $yyy, $zzz, $yyy, $zzz, $nero);
}

// calcolo accelerazioni

$vx = $v0["x"];
$vy = $v0["y"];
$vz = $v0["z"];

eval("\$ff = @($fx);");
$w0["x"] = $ff / $m;
eval("\$ff = @($fy);");
$w0["y"] = $ff / $m;
eval("\$ff = @($fz);");
$w0["z"] = $ff / $m;

// calcolo velocità dopo $dt

$v1["x"] = $v0["x"] + $w0["x"] * $dt;
$v1["y"] = $v0["y"] + $w0["y"] * $dt;
$v1["z"] = $v0["z"] + $w0["z"] * $dt;

// calcolo posizioni dopo $dt

$u1["x"] = $u0["x"] + $v0["x"] * $dt;
$u1["y"] = $u0["y"] + $v0["y"] * $dt;
$u1["z"] = $u0["z"] + $v0["z"] * $dt;

// reset

$u0["x"] = $u1["x"]; $u0["y"] = $u1["y"]; $u0["z"] = $u1["z"];
$v0["x"] = $v1["x"]; $v0["y"] = $v1["y"]; $v0["z"] = $v1["z"];

}

$xfin = $u1["x"] - $R0x;
$yfin = $u1["y"] - $R0y;
$zfin = $u1["z"] - $R0z;

$xxfin = $ex1 * $xfin + $ex2 * $yfin + $ex3 * $zfin;
$yyfin = $ey1 * $xfin + $ey2 * $yfin + $ey3 * $zfin;
$zzfin = $ez1 * $xfin + $ez2 * $yfin + $ez3 * $zfin;

echo "--- risultati dopo iterazioni :";

echo "Valori finali : x = " . $xxfin . " y = " . $yyfin . " z = " . $zzfin .
" R = " . sqrt($xxfin*$xxfin+$yyfin*$yyfin+$zzfin*$zzfin) . "";

$vxfin = $v1["x"];
$vyfin = $v1["y"];
$vzfin = $v1["z"];

$vxxfin = $ex1 * $vxfin + $ex2 * $vyfin + $ex3 * $vzfin;
$vyyfin = $ey1 * $vxfin + $ey2 * $vyfin + $ey3 * $vzfin;
$vzzfin = $ez1 * $vxfin + $ez2 * $vyfin + $ez3 * $vzfin;

echo "Valori finali : vx = " . $vxxfin . " vy = " . $vyyfin . " vz = " . $vzzfin .
" v = " . sqrt($vxxfin*$vxxfin+$vyyfin*$vyyfin+$vzzfin*$vzzfin) . "";

// disegna immagine

$nomefile = $REMOTE_ADDR;
$nomefile = "grafico" .str_replace(".", "-", $nomefile). ".png";

ImagePng ($im, $nomefile);

echo " ";
echo " ";
echo " ";
echo " ";

echo "<a href='$nomefile'>Visualizza il grafico</a>";

</script>

</body>

</html>

</html>