SPEED
GET_FUNC_VALUE.f90
Go to the documentation of this file.
1! Copyright (C) 2012 The SPEED FOUNDATION
2! Author: Ilario Mazzieri
3!
4! This file is part of SPEED.
5!
6! SPEED is free software; you can redistribute it and/or modify it
7! under the terms of the GNU Affero General Public License as
8! published by the Free Software Foundation, either version 3 of the
9! License, or (at your option) any later version.
10!
11! SPEED is distributed in the hope that it will be useful, but
12! WITHOUT ANY WARRANTY; without even the implied warranty of
13! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
14! Affero General Public License for more details.
15!
16! You should have received a copy of the GNU Affero General Public License
17! along with SPEED. If not, see <http://www.gnu.org/licenses/>.
18
33
34
35
36 real*8 function get_func_value(nb_fnc, type_fnc, ind_fnc, &
37 data_fnc, nb_data_fnc, id_fnc, time, dist,vel,timeit)
38
39
40
41 use binarysearch
42
43 implicit none
44
45 integer*4 :: nb_fnc,id_fnc,i,nb_data_fnc, idx !nb_timeval
46 integer*4 :: ind_start, ind_end, timeit
47
48 integer*4, dimension(nb_fnc) :: type_fnc
49 integer*4, dimension(nb_fnc +1) :: ind_fnc
50
51 real*8 :: pi,t_t0,t0,t1,v0,v1,omega, fp, fac
52 real*8 :: tau,scaling,hdur
53 real*8 :: amp, ps0, tplus, alpha,time,beta2,dist,vel
54
55 real*8, dimension(nb_data_fnc) :: data_fnc
56 real*8, dimension(1) :: valmax
57
58
59 get_func_value = 0.0d0
60
61 pi = 4.0d0 * datan(1.0d0)
62
63
64 ! write(*,*) type_fnc(id_fnc)
65 ! write(*,*) nb_fnc
66 ! write(*,*) type_fnc
67 ! write(*,*) ind_fnc
68 ! write(*,*) data_fnc
69 ! write(*,*) nb_data_fnc
70 ! write(*,*) id_fnc
71 ! read(*,*)
72
73 select case (type_fnc(id_fnc))
74
75 case(0)
76 get_func_value = 1.0d0
77
78 case(1)
79 t_t0 = time - data_fnc(ind_fnc(id_fnc) +1)
80 get_func_value = (1.0d0 - 2.0d0*data_fnc(ind_fnc(id_fnc))*t_t0*t_t0) &
81 * dexp(-1.0d0*data_fnc(ind_fnc(id_fnc))*t_t0*t_t0)
82
83 case(2)
84 pi = 4.0d0 * datan(1.0d0)
85 t_t0 = time - data_fnc(ind_fnc(id_fnc) +1)
86 get_func_value = dcos(pi*data_fnc(ind_fnc(id_fnc))*t_t0) &
87 * dexp(-0.5d0*data_fnc(ind_fnc(id_fnc)) &
88 * data_fnc(ind_fnc(id_fnc))*t_t0*t_t0)
89
90 case(3)
91 ind_start = ind_fnc(id_fnc); ind_end = ind_fnc(id_fnc+1)-3;
92
93! write(*,*) id_fnc, ind_fnc
94! read(*,*)
95! print*, ind_start, data_fnc(ind_start)
96! print*, ind_end, data_fnc(ind_end-1)
97 ! read*
98 ! print*, data_fnc(ind_start:ind_end:2)
99 ! read*
100
101 ! open(300,file='th.out',position='APPEND')
102 ! do i = ind_start, ind_end, 2
103 ! write(300,*) data_fnc(i), data_fnc(i+1)
104 ! enddo
105 ! close(300)
106
107
108 !valmax = maxval(data_fnc(ind_start:ind_end:2))
109 if (time >= data_fnc(ind_end-1)) then
110 v1 = data_fnc(ind_end + 2);
111 get_func_value = v1;
112 else
113 idx = binarysearch_real(data_fnc(ind_start:ind_end:2), time)
114
115! write(*,*) idx
116
117 t0 = data_fnc(ind_start-1 + 2*idx-1); t1 = data_fnc(ind_start-1 + 2*idx+1)
118 v0 = data_fnc(ind_start-1 + 2*idx); v1 = data_fnc(ind_start-1 + 2*idx+2)
119
120 get_func_value = (v1 - v0) / (t1 - t0) * (time - t0) + v0
121
122 endif
123
124
125
126! do i = ind_fnc(id_fnc), ind_fnc(id_fnc+1) -3,2
127! t0 = data_fnc(i); t1 = data_fnc(i +2)
128! v0 = data_fnc(i +1); v1 = data_fnc(i +3)
129! if ((time.ge.t0) .and. (time .le. t1)) then
130! write(*,*) time, t0, t1, v0, v1, GET_FUNC_VALUE
131 ! open(300,file='th.out',position='APPEND')
132 ! write(300,*) time, GET_FUNC_VALUE
133 ! close(300)
134 ! read*
135
136! GET_FUNC_VALUE = (v1 - v0) / (t1 - t0) * (time - t0) + v0
137! val2 = (v1 - v0) / (t1 - t0) * (time - t0) + v0
138! if (abs(val2-val1) .ne. 0) then
139! print*, time, val1, val2
140! read*
141! endif
142! return
143! endif
144! enddo
145
146 case(773) !!! time series given in file AH
147 get_func_value = data_fnc(ind_fnc(id_fnc) + timeit)
148
149 case(4)
150 t_t0 = time - data_fnc(ind_fnc(id_fnc) +1)
151 beta2 = data_fnc(ind_fnc(id_fnc))
152 get_func_value = 2.0d0*beta2*t_t0 &
153 * (-3.0d0 + 2.0d0*beta2*t_t0*t_t0) &
154 * dexp(-beta2*t_t0*t_t0)
155
156 case(5)
157 t_t0 = time - data_fnc(ind_fnc(id_fnc) +1)
158 fp = data_fnc(ind_fnc(id_fnc));
159 alpha = 2.d0*(pi*fp)**2
160 fac = 2*pi*fp*dsqrt(dexp(1.d0))
161 get_func_value = fac*t_t0*dexp(-alpha*t_t0**2);
162
163! write(*,*) time, data_fnc(ind_fnc(id_fnc) +1), fp, alpha, fac, GET_FUNC_VALUE
164! read(*,*)
165
166
167 case(6)
168 pi = 4.0d0 * datan(1.0d0)
169 t_t0 = time - data_fnc(ind_fnc(id_fnc) +1)
170 get_func_value = (2.d0*pi*data_fnc(ind_fnc(id_fnc))*t_t0) &
171 * dexp(-0.5d0*4.d0*pi*pi*data_fnc(ind_fnc(id_fnc)) &
172 * data_fnc(ind_fnc(id_fnc))*t_t0*t_t0)
173
174
175 case(7)
176 t_t0 = time - data_fnc(ind_fnc(id_fnc) +1)
177 get_func_value = (6*data_fnc(ind_fnc(id_fnc)) &
178 -24*data_fnc(ind_fnc(id_fnc))**2 * (t_t0**2) &
179 +8.0d0*data_fnc(ind_fnc(id_fnc))**3 * (t_t0**4)) &
180 * dexp(-1.0d0*data_fnc(ind_fnc(id_fnc))*t_t0*t_t0)
181
182
183 case(8)
184
185 omega = data_fnc(ind_fnc(id_fnc))*2*pi
186 get_func_value = dcos(omega*time)
187
188 case(9)
189
190 omega = data_fnc(ind_fnc(id_fnc))*2*pi
191 get_func_value = dsin(omega*time)
192
193
194
195 case(12)
196 !------------------------------------------------
197 ! 12 - sigmf(t,[a c]) = amp*(1/(1+exp(-a*(t-c))))
198 !
199 !
200 !
201 ! |
202 ! |............*************************......amp
203 ! | ** :
204 ! | a * :
205 ! | * :
206 ! | ** :
207 ! 0******----------------------------------> Time
208 ! | |
209 ! |____c___| |
210 ! | | |
211 ! n*(1/f)+t0
212 !
213 t_t0 = time - data_fnc(ind_fnc(id_fnc) +2)
214 get_func_value = data_fnc(ind_fnc(id_fnc)) &
215 * (1/(1+exp(-data_fnc(ind_fnc(id_fnc) +1)*(t_t0))))
216
217 case(13)
218 ! - GRENOBLE BENCHMARK
219 tau = data_fnc(ind_fnc(id_fnc)); scaling = data_fnc(ind_fnc(id_fnc) +1)
220 t0 = 2.0 * tau; hdur = tau/2.0; t_t0 = time - t0
221 get_func_value = 0.5d0 * ( 1.0d0 + erf(scaling*(t_t0)/hdur) )
222
223 case(14)
224 ! - SCEC BENCHMARK
225 tau = data_fnc(ind_fnc(id_fnc));
226 tau = tau*0.25d0;
227 amp = data_fnc(ind_fnc(id_fnc) +1)
228 t_t0 = time
229 if (t_t0 .lt. 0.0d0) then
230 get_func_value = 0.0d0
231 elseif (t_t0 .eq. 0.0d0) then
232 get_func_value = 0.5d0
233 else
234 get_func_value = amp * (1.0d0 - (1 + t_t0/tau)*exp(-t_t0/tau))
235 endif
236
237 case(15)
238 ! - EXPLOSION
239 ps0 = data_fnc(ind_fnc(id_fnc))
240 tplus = data_fnc(ind_fnc(id_fnc) +1); alpha = data_fnc(ind_fnc(id_fnc) +2)
241 t_t0 = time - data_fnc(ind_fnc(id_fnc) +3)
242 if (t_t0 .lt. 0.0d0) then
243 get_func_value = 0.0d0
244 else
245 get_func_value = ps0 * (1 - (1 + time/tplus)*exp(-alpha*time/tplus))
246 endif
247
248! case(30)
249! t_t0 = time - dist/vel;
250! if (t_t0 .le. 0) then
251! GET_FUNC_VALUE = 0.d0;
252! else
253! GET_FUNC_VALUE = 2.0d0*beta2*t_t0 &
254! * (-3.0d0 + 2.0d0*beta2*t_t0*t_t0) &
255! * dexp(-beta2*t_t0*t_t0)
256! endif
257
258 case(30)
259
260
261 do i = ind_fnc(id_fnc),ind_fnc(id_fnc+1) -3,2
262 t0 = data_fnc(i) ; t1 = data_fnc(i +2)
263 v0 = data_fnc(i +1); v1 = data_fnc(i +3)
264
265 ! write(*,*) t0,t1,v0,v1
266 if (((time -dist/vel) .ge. t0) .and. ((time -dist/vel) .le. t1)) then
267 get_func_value = (v1 - v0) / (t1 - t0) * (time - t0) + v0
268 return
269 endif
270 enddo
271
272 ! ind_start = ind_fnc(id_fnc); ind_end = ind_fnc(id_fnc+1)-3;
273
274
275 ! if (time - dist/vel >= data_fnc(ind_end-1)) then
276 ! v1 = data_fnc(ind_end + 2);
277 ! GET_FUNC_VALUE = v1;
278 ! else
279 ! idx = binarySearch_real(data_fnc(ind_start:ind_end:2), time-dist/vel)
280 ! t0 = data_fnc(2*idx-1); t1 = data_fnc(2*idx+1)
281 ! v0 = data_fnc(2*idx); v1 = data_fnc(2*idx+2)
282 !
283 ! GET_FUNC_VALUE = (v1 - v0) / (t1 - t0) * (time - t0) + v0!
284 !
285 ! endif
286
287
288
289
290
291
292 case(60,62)
293 ! FUNCTION FOR G/G0
294 do i = ind_fnc(id_fnc),ind_fnc(id_fnc+1) -3,2
295
296 t0 = data_fnc(i); t1 = data_fnc(i +2)
297 v0 = data_fnc(i +1); v1 = data_fnc(i +3)
298
299 if (abs(time).le.data_fnc(ind_fnc(id_fnc))) then
300 get_func_value = data_fnc(ind_fnc(id_fnc)+1)
301 elseif ((abs(time).ge.t0).and.(abs(time).le.t1)) then
302 get_func_value = (v1 - v0) / (t1 - t0) * (abs(time) - t0) + v0
303 elseif (abs(time).ge.data_fnc(ind_fnc(id_fnc+1)-2)) then
304 get_func_value = data_fnc(ind_fnc(id_fnc+1)-1)
305 endif
306 enddo
307
308 case(61,63)
309 ! FUNCTION FOR DAMPING
310 do i = ind_fnc(id_fnc),ind_fnc(id_fnc+1) -3,2
311 t0 = data_fnc(i); t1 = data_fnc(i +2)
312 v0 = data_fnc(i +1); v1 = data_fnc(i +3)
313
314 if (abs(time).le.data_fnc(ind_fnc(id_fnc))) then
315 get_func_value = data_fnc(ind_fnc(id_fnc)+1)
316 elseif ((abs(time).ge.t0).and.(abs(time).le.t1)) then
317 get_func_value = (v1 - v0) / (t1 - t0) * (abs(time) - t0) + v0
318 elseif (abs(time).ge.data_fnc(ind_fnc(id_fnc+1)-2)) then
319 get_func_value = data_fnc(ind_fnc(id_fnc+1)-1)
320 endif
321 enddo
322
323 case(99)
324 ! - CASHIMA1 BENCHMARK
325 tau = data_fnc(ind_fnc(id_fnc)); amp = data_fnc(ind_fnc(id_fnc) +1)
326 t_t0 = time
327 get_func_value = ( exp( - (((2.0d0 * 2.0d0*dasin(1.0d0))*1.5d0)&
328 *(t_t0 - tau)/amp)**2.0d0)&
329 * cos(((2.0d0 * 2.0d0*dasin(1.0d0))*1.5d0)&
330 *(t_t0 - tau) + 2.0d0*dasin(1.0d0)/2.0d0))
331
332 case(100)
333 ! - TESTMODE
334 amp = data_fnc(ind_fnc(id_fnc))
335 get_func_value = dsin(amp*pi*time)
336
337
338 case(101)
339 get_func_value = time
340 case default
341 get_func_value = 0.d0
342
343 end select
344
345
346 return
347
348 end function get_func_value
349
get_func_value
real *8 function get_func_value(nb_fnc, type_fnc, ind_fnc, data_fnc, nb_data_fnc, id_fnc, time
Computes time evolution function.
Definition GET_FUNC_VALUE.f90:38
binarysearch::binarysearch_real
pure recursive integer function, public binarysearch_real(vec, scal, min, max)
Definition MODULES.f90:679