Programma Fortran per equazione moto, campo elettromagnetico

[alexsia1]

Buonasera a tutti, studio Fisica all'università ed è la prima volta che scrivo nel vostro forum, dunque spero di esprimere al meglio le mie perplessità.

Dovrei creare un programma in linguaggio Fortran che descriva il moto di una carica se applichiamo su di essa la forza di Lorentz F=q*(E+vxB), (con x sto indicando il prodotto vettoriale).
Ho creato un programma con l'algoritmo di Eulero (così come richiesto dal problema) che viene compilato correttamente, ma fornisce dei dati in output che, se graficati con GnuPlot, danno sempre una retta. Succede anche se modifico i dati e la cosa mi rende piuttosto perplessa.
Dove sto sbagliando?

Qui di seguito vi scrivo il testo del problema e il mio programma, in modo da rendervi più chiaro il discorso.

Testo del problema:
Usare l' algoritmo di Eulero per integrare le equazioni del moto di una
particella carica in campi elettrici e magnetici statici.

La legge di forza e'

F = q ( E + v x B )

Dove E e B sono dei vettori (campo elettrico e magnetico) e "x"
rappresenta il prodotto vettoriale. L' esercizio non dipende dalla
conoscenza dell' elettromagnetismo, trattandosi di campi assegnati; ci
chiediamo solo quale sara' la conseguenza sul moto di un punto materiale
della legge di forza sopra data. L' unica informazione supplementare utile
per lo svolgimento dell' esercizio puo' essere che la componente magnetica
della forza (q v x B ) essendo sempre ortogonale allo spostamento non
compie lavoro e quindi l' energia totale si conserva. L' energia
potenziale U, che e' dovuta al solo termine col campo elettrico, e'
ottenibile sapendo
che la componente i-esima della forza e' data da - dU/dx_i, x_i essendo
la i-esima coordinata cartesiana).

Discutere la scelta dello step temporale di integrazione e verificare il
comportamento del programma nei seguenti casi facendo partire dall' origine
al tempo t=0 una particella di carica e massa unitaria ( indichero' le
componenti cartesiane x,y,z di un vettore mediante terne
ordinate):

caso 1)

E = ( 0,0,0)

B = ( 0,0,1)

v(t=0) = (0,1,1)

caso 2)

E = ( 0,0,1)

B = ( 0,0,1)

v(t=0) = (0,1,0)

caso 3)

E = ( cos(2*pi*x),0,0)

B = ( 1,0,0)

v(t=0) = (0,1,0)

Soluzione proposta:

program j_san								! Programma per lo studio del moto di una carica q;
									! come conseguenza della Forza di Lorentz;
	implicit none
! Definiamo i tipi dati del problema
	integer :: i, nstep, error
	real, PARAMETER :: m=1.0, q=1.0
	real :: t=0.0, dt, dx, dv, pi					! emec variabile facoltativa, inserita solo per maggiore chiarezza;
	real, DIMENSION(3) :: F, B, pos, v, vxB, x, y, ekin, epot, emec	! vxB rappresenta il prodotto vettoriale.
	real, ALLOCATABLE, DIMENSION(:) :: E

! Definiamo le variabili
	pos = (/0.0,0.0,0.0/)
	print*,"Imposta l'intervallo d'integrazione dt più adeguato per questo problema:"
	read*,dt
	print*,"Imposta il numero di volte che vuoi ripetere il ciclo:"
	read*,nstep
!Primo caso:
	print*,"Moto della particella nel primo caso"
	ALLOCATE(E(3),STAT=error)
	do i = 1,nstep
		t = t + dt
		E = (/0.0,0.0,0.0/); B = (/0.0,0.0,1.0/); v = (/0.0,1.0,1.0/)
		vxB = (/(v(2)*B(3) - v(3)*B(2)),(v(3)*B(1) - v(1)*B(3)),(v(1)*B(2) - v(2)*B(1))/)
		F = q*(E + vxB)
		pos = pos + v*t + 0.5 * (F/m) * t**2
		v = v + (F/m)*t
		x = v*t
		y = 0.5 * q/m * E * t**2
		epot = - F*x
		ekin = 0.5 * m * v**2
	   	emec = ekin + epot
		print*,"i=",i,"t=",t,"pos=(",pos,")","v=(",v,")","x=(",x,")","y=(",y,")","ekin=(",ekin,")","epot=(",epot,")"
	end do
!Secondo caso:
	print*,"Moto della particella nel secondo caso"
	ALLOCATE(E(3),STAT=error)
	do i = 1,nstep
		t = t + dt
		E = (/0.0,0.0,1.0/); B = (/0.0,0.0,1.0/); v = (/0.0,1.0,0.0/)
		vxB = (/(v(2)*B(3) - v(3)*B(2)),(v(3)*B(1) - v(1)*B(3)),(v(1)*B(2) - v(2)*B(1))/)
		F = q*(E + vxB)
		pos = pos + v*t + 0.5 * (F/m) * t**2
		v = v + (F/m)*t
		x = v*t
		y = 0.5 * q/m * E * t**2
		epot = - F*x
		ekin = 0.5 * m * v**2
		emec = ekin + epot
		print*,"i=",i,"t=",t,"pos=(",pos,")","v=(",v,")","x=(",x,")","y=(",y,")","ekin=(",ekin,")","epot=(",epot,")"
	end do
!Terzo caso:
	print*,"Moto della particella nel terzo caso"
	ALLOCATE(E(5),STAT=error)
	do i = 1,nstep
		t = t + dt
		E = (/cos(2*pi*x),0.0,0.0/); B = (/1.0,0.0,0.0/); v = (/0.0,1.0,0.0/)
		vxB = (/(v(2)*B(3) - v(3)*B(2)),(v(3)*B(1) - v(1)*B(3)),(v(1)*B(2) - v(2)*B(1))/)
		F = q*(E + vxB)
		pos = pos + v*t + 0.5 * (F/m) * t**2
		v = v + (F/m)*t
		x = v*t
		y = 0.5 * q/m * E * t**2
		epot = - F*x
		ekin = 0.5 * m * v**2
		emec = ekin + epot
		print*,"i=",i,"t=",t,"pos=(",pos,")","v=(",v,")","x=(",x,")","y=(",y,")","ekin=(",ekin,")","epot=(",epot,")"
	end do

end program j_san

Mi scuso in anticipo per eventuali errori che per voi saranno magari inaccettabili, ma mi sono appena affacciata nel mondo dei linguaggi di programmazione.
Grazie in anticipo a chi mi sarà di aiuto.

[apatriarca]

Ciao alexsia,

benvenuta nel forum. Mi sono preso la libertà di editare il tuo post in modo da mostrare il codice formattato.

Venendo al codice.. Non mi è chiaro il significato delle variabili vettoriali \(x\) e \(y\). La velocità viene sovrascritta con il valore iniziale tutte le iterazioni invece di inizializzarla al di fuori del ciclo e poi lasciare che vari in base alla forza che hai calcolato nella tua iterazione. La \(x\) nel terzo caso credo sia semplicemente la seconda componente della posizione. Ma la tua \(x\) è calcolata diversamente.

[alexsia1]

Ciao, grazie per aver editato al posto mio
Quindi consigli di eliminare le variabili x e y? In realtà le avevo introdotte perché volevo imporre l'equazione della parabola (un possibile caso di equazione del moto in campo elettrico) e vedere cosa succedeva, in un modo un po' grezzo e confusionario, me ne rendo conto solo adesso.
Ora provo a modificare la velocità e vedo cosa mi esce.
Grazie!

[alexsia1]

Ho apportato i cambiamenti da te consigliati, adesso nei primi due casi ho una parabola e nel terzo una retta. Come mai fornisce questi risultati?
Scusa se chiedo ogni dettaglio, ma sono abbastanza confusa..

[apatriarca]

La parabola non è l'unica soluzione possibile al tuo problema. Anche una retta è una soluzione possibile se il campo magnetico è parallelo alla velocità lungo tutto il moto. Ed è in effetti questo il caso nel terzo punto. Hai infatti che il campo magnetico \(B\) è uguale a zero su tutto il piano \(x=0\) nel quale ti trovi nell'istante \(t=0\). Siccome la componente della velocità nella direzione \(x\) è nulla, rimarrai sempre in tale piano muovendoti lungo la retta \(y = t\).

[alexsia1]

Perfetto! Ti ringrazio tantissimo!