-
Notifications
You must be signed in to change notification settings - Fork 0
/
fractal.f90
179 lines (137 loc) · 6.39 KB
/
fractal.f90
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
subroutine fractal
!==============================================================!
! Esta subrutina guada los datos del fractal. !
!==============================================================!
!===============================================================
! Usamos el módulo 'vars' para tener información de los parámetros de entrada.
use vars
! Por precaución no dejamos ningún tipo implícito.
implicit none
! Declaramos variables útiles.
integer i, j, k
real(8) tg, tl
real(8), allocatable, dimension (:) :: tiempo_g, tiempo_l
allocate( tiempo_g(Nx*Nx), tiempo_l(Nx*Nx) )
!===============================================================
! Definimos las condiciones iniciales.
Th1(1) = -pi
Th2(1) = -pi
Xi1(1) = cero
Xi2(1) = cero
t(1) = cero
!===============================================================
! Abrimos un archivo para guardar los datos, accedemos a él con la etiqueta '10'.
open(20, file = 'fractal.csv', status = 'replace')
! Escribimos el cabezal del archivo.
write(20,*) '#Th1, Th2, tg, tl'
! Imprimimos en terminal información útil.
write(*,*) '---------------------------------'
write(*,*) '| Vamos a comenzar a solucionar |'
write(*,*) '---------------------------------'
write(*,*) '| Línea | Total |'
write(*,*) '---------------------------------'
! Iniciamos el contador de líneas y el contador de índices para los tiempos.
j = 1
k = 0
! Llenamos el arreglo temporal.
do i=1, Nt-1
t(i+1) = t(i) + dt
end do
!===============================================================
! Iniciamos los ciclos.
do while (Th1(1) <= pi)
do while (Th2(1) <= pi)
!---------------------------------------------------------------
! Solucionamos solo la mitad izquierda.
if (Th1(1) <= 0) then
! Esta es la condición de giro.
if ( alpha*beta*cos(Th2(1)) + (1+alpha)*cos(Th1(1)) > 1 + alpha*(1-beta) ) then
!---------------------------------------------------------------
! Si no gira ni latigea podemos decir que lo hace al tiempo final.
tg = tf
tl = tf
!---------------------------------------------------------------
! Escribimos la información en el archivo.
write(20,200) Th1(1), Th2(1), tg, tl
!---------------------------------------------------------------
! Guardamos el tiempo de giro y latiganzo en un arreglo.
tiempo_g(k) = tg
tiempo_l(k) = tl
else
!---------------------------------------------------------------
! Si no se cumple la condición calculemos en qué momento gira.
! Hacemos el método de Runge-Kutta en sí.
do i=1, Nt-1
!-------------------------------------------
call rk4(Th1(i), Th2(i), Xi1(i), Xi2(i))
!-------------------------------------------
Th1(i+1) = Th1(i) + ( k1Th1 + dos*k2Th1 + dos*k3Th1 + k4Th1 )*( dt*sexto )
Th2(i+1) = Th2(i) + ( k1Th2 + dos*k2Th2 + dos*k3Th2 + k4Th2 )*( dt*sexto )
Xi1(i+1) = Xi1(i) + ( k1Xi1 + dos*k2Xi1 + dos*k3Xi1 + k4Xi1 )*( dt*sexto )
Xi2(i+1) = Xi2(i) + ( k1Xi2 + dos*k2Xi2 + dos*k3Xi2 + k4Xi2 )*( dt*sexto )
end do
!---------------------------------------------------------------
! Reiniciamos los tiempos de giro y latiganzo.
tg = tf
tl = tf
!---------------------------------------------------------------
! Calculamos el tiempo de giro.
do i = 1, Nt
if ( ( abs(Th1(i)) >= pi ) .or. ( abs(Th2(i)) >= pi ) ) then
tg = t(i-1)
exit
end if
end do
!---------------------------------------------------------------
! Calculamos el tiempo de latiganzo.
do i = 1, Nt
if ( Xi1(i)*Xi2(i) <= 0 ) then
tl = t(i-1)
exit
end if
end do
!---------------------------------------------------------------
! Escribimos la información en el archivo.
write(20,200) Th1(1), Th2(1), tg, tl
!---------------------------------------------------------------
! Guardamos el tiempo de giro y latiganzo en un arreglo.
tiempo_g(k) = tg
tiempo_l(k) = tl
end if
!---------------------------------------------------------------
! Actualizamos el contador de los tiempos.
k = k + 1
!---------------------------------------------------------------
! Hacemos la reflexión para la mitad derecha.
else
!---------------------------------------------------------------
! Calculamos los tiempos de giro y latiganzo en función de los ya guardados.
tg = tiempo_g(k)
tl = tiempo_l(k)
!---------------------------------------------------------------
! Escribimos en el archivo.
write(20,200) Th1(1), Th2(1), tg, tl
!---------------------------------------------------------------
! Movemos el contador de los tiempos al contrario.
k = k - 1
end if
! Damos un paso hacia arriba en la línea para \theta_2(0).
Th2(1) = Th2(1) + dx
end do
! Damos un paso en la malla para \theta_1(0) y reiniciamos el valor de \theta_2(0).
Th1(1) = Th1(1) + dx
Th2(1) = -pi
! Imprimimos en pantalla en qué línea vamos.
write(*,"(a5,i7,a9,i7,a6)") ' | ',j,' de ',Nx,' |'
! Actualizamos el contador de líneas.
j = j + 1
end do
! Imprimimos información útil en pantalla.
write(*,*) '---------------------------------'
write(*,*) '| Ya terminamos, Camarada UwUr |'
write(*,*) '---------------------------------'
! Cerramos el archivo donde escribimos la información.
close(20)
! Definimos el formato de escritura en el archivo, este es estilo cvs.
200 format (f20.16,',',f20.16,',',f20.16,',',f20.16)
end subroutine fractal