COSMOS v7.655  COSMOSv7655
(AirShowerMC)
timeslice.f File Reference

Go to the source code of this file.

Functions/Subroutines

program __timeslice.f__
 
subroutine slice (mint, maxt, dt)
 
subroutine thinning (i1, i2)
 
subroutine mvdata (f, j, nc)
 
subroutine mvdatai (f, j, nc)
 
subroutine mkfiles
 
subroutine mkskelfiles
 
subroutine mkafile
 
subroutine mkaskelfile
 
subroutine wtaskel (fno, xa, ya, za, code, chg, np, mulalpha)
 
subroutine code2rgb (codei, chgi, r, g, b, alpha)
 

Function/Subroutine Documentation

◆ __timeslice.f__()

program __timeslice.f__ ( )

Definition at line 4 of file timeslice.f.

References copenfw2(), klena(), mkafile(), mkaskelfile(), mkfiles(), mkskelfiles(), n, slice(), and thinning().

4  character*120 tracefile
Here is the call graph for this function:

◆ code2rgb()

subroutine code2rgb ( integer*2  codei,
integer*2  chgi,
real  r,
real  g,
real  b,
real  alpha 
)

Definition at line 381 of file timeslice.f.

References code, parameter(), and true.

Referenced by wtaskel().

381  implicit none
382  integer*2 codei, chgi
383  real r, g, b, alpha
384 
385  integer ncolor, mncolor
386  parameter( mncolor = 17 )
387  type colortab
388  integer*2 code
389  integer*2 chg
390  real r
391  real g
392  real b
393  real alpha
394  end type colortab
395  type(colortab):: tab(mncolor)
396 
397  character*80 input
398  integer i
399 
400  integer first
401  data first/0/
402  save first, ncolor
403 
404 
405  if(first .eq. 0) then
406  input = ' '
407  open(13, file='colortab')
408  i = 0
409  do while(.true.)
410  read(13,'(a)',end=10) input
411  if(input(1:1) .ne. "#" .and. input .ne. ' ') then
412  i = i + 1
413  if(i .gt. mncolor) then
414  write(0,*) ' too many color spec. in colortab'
415  stop 99999
416  endif
417  read(input,*)
418  * tab(i).code, tab(i).chg,
419  * tab(i).r, tab(i).g, tab(i).b,
420  * tab(i).alpha
421  endif
422  enddo
423  10 continue
424  ncolor = i
425  first = 1
426  close(13)
427  endif
428 
429  do i = 1, ncolor
430  if(tab(i).code .eq. codei .and. tab(i).chg .eq. chgi) goto 100
431  enddo
432  i = ncolor
433  100 continue
434  r = tab(i).r
435  g = tab(i).g
436  b = tab(i).b
437  alpha = tab(i).alpha
integer npitbl real *nx parameter(n=101, npitbl=46, nx=n-1) real *8 uconst
*Zfirst p fm *Zfirst p Zfirst p Zfirst p *Zfirst p *Zfirst pos xyz r
Definition: ZavoidUnionMap.h:1
logical, save first
Definition: cNRLAtmos.f:8
nodes i
block data cblkElemag data *AnihiE ! Eposi< 1 TeV, anihilation considered *X0/365.667/, ! radiation length of air in kg/m2 *Ecrit/81.e-3/, ! critical energy of air in GeV *MaxComptonE/1./, ! compton is considered below 1 GeV *MaxPhotoE/1.e-3/, ! above this, PhotoElectric effect neg. *MinPhotoProdE/153.e-3/, ! below 153 MeV, no gp --> hadrons ! scattering const not MeV *Knockon true
Definition: cblkElemag.h:7
*Zfirst p fm *Zfirst p Zfirst p code
Definition: ZavoidUnionMap.h:1
dE dx *! Nuc Int sampling table g
Definition: cblkMuInt.h:130
real(4), save b
Definition: cNRLAtmos.f:21
Here is the call graph for this function:
Here is the caller graph for this function:

◆ mkafile()

subroutine mkafile ( )

Definition at line 301 of file timeslice.f.

References idx(), x, y, and z.

Referenced by __timeslice.f__(), __timeslice2.f__(), __timesliceg77.f__(), and slice1tbin().

301  implicit none
302  include "Zprivate2.f"
303 
304  integer i, j
305  character* 130 filename
306 
307  integer filec
308  integer klena
309  data filec/0/
310  save filec
311 
312  if( filec .eq. 0 ) then
313  filename = dir(1:klena(dir))//"/timesorted.dat"
314  open(20, file=filename,
315  * form='formatted')
316  filec = 1
317  endif
318 
319  do j = jmin, jmax
320  do i= 1, idx(j)
321  write(20,'(3g16.8, i3,i3)' )
322  * x(i, j), y(i,j), z(i,j), codex(i,j), chgx(i,j)
323  enddo
324  if( idx(j) .gt. 0 ) write(20,*)
325  enddo
nodes z
nodes i
averg real MaxCPU integer idx(Maxp)
********************block data cblkHeavy ********************integer j data *HeavyG2symbol *data *HeavyG2code kiron data j
Definition: cblkHeavy.h:36
latitude latitude this system is used *****************************************************************! type coord sequence union map real y
Definition: Zcoord.h:25
integer function klena(cha)
Definition: klena.f:20
! structure defining a particle at production ! Basic idea of what is to be contained in ! the particle structue is that dynamical ones should be included those derivable from the particle code ! is not included ******************************************************type fmom momentum sequence union map real e endmap map real * x
Definition: Zptcl.h:21
Here is the call graph for this function:
Here is the caller graph for this function:

◆ mkaskelfile()

subroutine mkaskelfile ( )

Definition at line 328 of file timeslice.f.

References idx(), wtaskel(), x, y, and z.

Referenced by __timeslice.f__(), __timeslice2.f__(), __timesliceg77.f__(), and slice1tbin().

328  implicit none
329  include "Zprivate2.f"
330 
331  integer i, j
332  character* 130 filename
333 
334  integer filec
335  integer klena
336  data filec/0/
337  save filec
338 
339  if( filec .eq. 0 ) then
340  filename = dir(1:klena(dir))//"/timesorted.skel"
341  open(20, file=filename,
342  * form='formatted')
343  filec = 1
344  endif
345 
346  do j = jmin, jmax
347  if(idx(j) .gt. 0) then
348  write(20,'(a)') "SKEL"
349  call wtaskel(20,
350  * x(1, j), y(1,j), z(1,j), codex(1,j), chgx(1,j), idx(j),
351  * mulalpha)
352  write(20,*)
353  endif
354  enddo
nodes z
nodes i
subroutine wtaskel(fno, xa, ya, za, code, chg, np, mulalpha)
Definition: timeslice.f:359
averg real MaxCPU integer idx(Maxp)
********************block data cblkHeavy ********************integer j data *HeavyG2symbol *data *HeavyG2code kiron data j
Definition: cblkHeavy.h:36
latitude latitude this system is used *****************************************************************! type coord sequence union map real y
Definition: Zcoord.h:25
integer function klena(cha)
Definition: klena.f:20
! structure defining a particle at production ! Basic idea of what is to be contained in ! the particle structue is that dynamical ones should be included those derivable from the particle code ! is not included ******************************************************type fmom momentum sequence union map real e endmap map real * x
Definition: Zptcl.h:21
Here is the call graph for this function:
Here is the caller graph for this function:

◆ mkfiles()

subroutine mkfiles ( )

Definition at line 242 of file timeslice.f.

References idx(), x, y, and z.

Referenced by __timeslice.f__(), __timeslice2.f__(), __timesliceg77.f__(), and slice1tbin().

242  implicit none
243  include "Zprivate2.f"
244 
245  character* 130 filename
246 
247  integer filec
248  integer klena
249  data filec/0/
250  save filec
251 
252  integer i, j
253 
254  do j = jmin, jmax
255  if(idx(j) .gt. 0) then
256  filec = filec + 1
257  write(filename,'(a, a, i5.5,a)') dir(1:klena(dir)),
258  * "/ts", j, ".dat"
259 ! * "/ts", filec, ".dat"
260  open(20, file=filename, form='formatted')
261  do i= 1, idx(j)
262  write(20,'(3g16.8,i3,i3)')
263  * x(i, j), y(i,j), z(i,j), codex(i,j), chgx(i,j)
264  enddo
265  close(20)
266  endif
267  enddo
nodes z
nodes i
averg real MaxCPU integer idx(Maxp)
********************block data cblkHeavy ********************integer j data *HeavyG2symbol *data *HeavyG2code kiron data j
Definition: cblkHeavy.h:36
latitude latitude this system is used *****************************************************************! type coord sequence union map real y
Definition: Zcoord.h:25
integer function klena(cha)
Definition: klena.f:20
! structure defining a particle at production ! Basic idea of what is to be contained in ! the particle structue is that dynamical ones should be included those derivable from the particle code ! is not included ******************************************************type fmom momentum sequence union map real e endmap map real * x
Definition: Zptcl.h:21
Here is the call graph for this function:
Here is the caller graph for this function:

◆ mkskelfiles()

subroutine mkskelfiles ( )

Definition at line 271 of file timeslice.f.

References idx(), wtaskel(), x, y, and z.

Referenced by __timeslice.f__(), __timeslice2.f__(), __timesliceg77.f__(), and slice1tbin().

271  implicit none
272  include "Zprivate2.f"
273 
274  character* 130 filename
275 
276  integer filec
277  integer klena
278  data filec/0/
279  save filec
280 
281  integer i, j
282 
283  do j = jmin, jmax
284  if(idx(j) .gt. 0) then
285  filec = filec + 1
286  write(filename,'(a, a, i5.5, a)') dir(1:klena(dir)),
287  * "/ts", j, ".skel"
288 ! * "/ts", filec, ".skel"
289  open(20, file=filename, form='formatted')
290 
291  write(20, '(a)') "SKEL"
292  call wtaskel(20, x(1,j), y(1,j), z(1,j),
293  * codex(1,j), chgx(1, j), idx(j), mulalpha)
294 
295  close(20)
296  endif
297  enddo
nodes z
nodes i
subroutine wtaskel(fno, xa, ya, za, code, chg, np, mulalpha)
Definition: timeslice.f:359
averg real MaxCPU integer idx(Maxp)
********************block data cblkHeavy ********************integer j data *HeavyG2symbol *data *HeavyG2code kiron data j
Definition: cblkHeavy.h:36
latitude latitude this system is used *****************************************************************! type coord sequence union map real y
Definition: Zcoord.h:25
integer function klena(cha)
Definition: klena.f:20
! structure defining a particle at production ! Basic idea of what is to be contained in ! the particle structue is that dynamical ones should be included those derivable from the particle code ! is not included ******************************************************type fmom momentum sequence union map real e endmap map real * x
Definition: Zptcl.h:21
Here is the call graph for this function:
Here is the caller graph for this function:

◆ mvdata()

subroutine mvdata ( real*4, dimension(maxp, n f,
integer  j,
integer  nc 
)

Definition at line 191 of file timeslice.f.

References idx().

Referenced by thinning().

191  implicit none
192  include "Zprivate2.f"
193 
194  integer j
195  real*4 f(maxp, n)
196 
197  real*4 temp(maxp)
198  integer i, nc, l
199 
200  nc = 0
201 
202  do i = 1, idx(j)
203  l = si(i)
204  if(l .gt. 0 ) then
205  nc = nc + 1
206  temp(nc) = f(l, j)
207  endif
208  enddo
209 
210  do i = 1, nc
211  f(i, j) = temp(i)
212  enddo
213 
214 
nodes i
averg real MaxCPU integer idx(Maxp)
********************block data cblkHeavy ********************integer j data *HeavyG2symbol *data *HeavyG2code kiron data j
Definition: cblkHeavy.h:36
real(4), dimension(:), allocatable, save temp
Definition: cNRLAtmos.f:29
!onst int maxp
Definition: Zprivate.h:3
real cut integer nc
Definition: Zprivate.h:1
integer n
Definition: Zcinippxc.h:1
dE dx *! Nuc Int sampling table f
Definition: cblkMuInt.h:130
Here is the call graph for this function:
Here is the caller graph for this function:

◆ mvdatai()

subroutine mvdatai ( integer*2, dimension(maxp, n f,
integer  j,
integer  nc 
)

Definition at line 217 of file timeslice.f.

References idx().

Referenced by thinning().

217  implicit none
218  include "Zprivate2.f"
219  integer j
220  integer*2 f(maxp, n)
221 
222  integer*2 temp(maxp)
223  integer i, nc, l
224  nc = 0
225 
226  do i = 1, idx(j)
227  l = si(i)
228  if(l .gt. 0 ) then
229  nc = nc + 1
230  temp(nc) = f(l, j)
231  endif
232  enddo
233 
234  do i = 1, nc
235  f(i, j) = temp(i)
236  enddo
237 
238 
nodes i
averg real MaxCPU integer idx(Maxp)
********************block data cblkHeavy ********************integer j data *HeavyG2symbol *data *HeavyG2code kiron data j
Definition: cblkHeavy.h:36
real(4), dimension(:), allocatable, save temp
Definition: cNRLAtmos.f:29
!onst int maxp
Definition: Zprivate.h:3
real cut integer nc
Definition: Zprivate.h:1
integer n
Definition: Zcinippxc.h:1
dE dx *! Nuc Int sampling table f
Definition: cblkMuInt.h:130
Here is the call graph for this function:
Here is the caller graph for this function:

◆ slice()

subroutine slice ( real*4  mint,
real*4  maxt,
real*4  dt 
)

Definition at line 69 of file timeslice.f.

References idx(), maxp, n, thinning(), true, x, y, and z.

Referenced by __timeslice.f__(), __timeslice2.f__(), __timesliceg77.f__(), and slice1tbin().

69  implicit none
70  include "Zprivate2.f"
71 
72  real*4 mint, maxt, dt
73  integer i, j
74  real*4 t, k
75  integer*2 code, chg
76  integer it1, it2
77  real*4 xx1, yy1, zz1, time1
78  real*4 xx2, yy2, zz2, time2
79 
80  character* 120 input
81 
82  do j = 1, n
83  idx(j) = 0
84  enddo
85 
86  do while (.true.)
87  input = ' '
88  read(11, '(a)', end=1000 ) input
89  if( input .ne. " " ) then
90  read(input, * ) xx1, yy1, zz1, time1, code, chg
91  do while( input .ne. " ")
92  read(11, '(a)', end =1000) input
93  if( input .ne. " " ) then
94  read(input, * ) xx2, yy2, zz2, time2, code, chg
95  it1 = (time1- mint)/dt +1
96  it2 = (time2- mint)/dt
97  if(it1 .lt. 1) goto 900
98  if(it2 .gt. n) goto 900
99  if(time2 .gt. maxt) goto 900
100  jmin = min(jmin, it1)
101  jmax = max(jmax, it2)
102  do j = it1, it2
103  t = j*dt + mint
104  if( idx(j) .ge. maxp ) then
105  write(0,*)
106  * 'ptcls at time=', t , ' > ', maxp
107  write(0,*)
108  * 'try to thin ptcls in the same pixel'
109  call thinning(j,j)
110  endif
111  idx(j) = idx(j) + 1
112  if(idx(j) .gt. maxp) then
113  write(0,*) 'ptcls > ', maxp
114  write(0,*) 'thinning not worked'
115  stop 1234
116  endif
117 
118 ! (x2-x1)/(t2-t1) *(t-t1) +x1
119  k = (t - time1) /(time2-time1)
120  x(idx(j), j) = (xx2-xx1) * k + xx1
121  y(idx(j), j) = (yy2-yy1) * k + yy1
122  z(idx(j), j) = (zz2-zz1) * k + zz1
123  codex(idx(j), j) = code
124  chgx(idx(j), j) = chg
125  enddo
126  900 continue
127  xx1 = xx2
128  yy1 = yy2
129  zz1 = zz2
130  time1 = time2
131  endif
132  enddo
133  endif
134  enddo
135  1000 continue
nodes z
nodes i
subroutine thinning(i1, i2)
Definition: timeslice.f:141
block data cblkElemag data *AnihiE ! Eposi< 1 TeV, anihilation considered *X0/365.667/, ! radiation length of air in kg/m2 *Ecrit/81.e-3/, ! critical energy of air in GeV *MaxComptonE/1./, ! compton is considered below 1 GeV *MaxPhotoE/1.e-3/, ! above this, PhotoElectric effect neg. *MinPhotoProdE/153.e-3/, ! below 153 MeV, no gp --> hadrons ! scattering const not MeV *Knockon true
Definition: cblkElemag.h:7
averg real MaxCPU integer idx(Maxp)
********************block data cblkHeavy ********************integer j data *HeavyG2symbol *data *HeavyG2code kiron data j
Definition: cblkHeavy.h:36
*Zfirst p fm *Zfirst p Zfirst p code
Definition: ZavoidUnionMap.h:1
latitude latitude this system is used *****************************************************************! type coord sequence union map real y
Definition: Zcoord.h:25
nodes t
!onst int maxp
Definition: Zprivate.h:3
integer n
Definition: Zcinippxc.h:1
! structure defining a particle at production ! Basic idea of what is to be contained in ! the particle structue is that dynamical ones should be included those derivable from the particle code ! is not included ******************************************************type fmom momentum sequence union map real e endmap map real * x
Definition: Zptcl.h:21
Here is the call graph for this function:
Here is the caller graph for this function:

◆ thinning()

subroutine thinning ( integer  i1,
integer  i2 
)

Definition at line 141 of file timeslice.f.

References false, idx(), kqsortr(), mvdata(), mvdatai(), true, x, y, and z.

Referenced by __timeslice.f__(), __timeslice2.f__(), __timesliceg77.f__(), slice(), and slice1tbin().

141  implicit none
142  include "Zprivate2.f"
143 ! culling overlapping points in a pixel
144 ! this is very rough thinning
145 ! sorr each x
146  integer i1, i2
147  integer nc, j, i, k1, k2, l, nc0
148  logical dothin
149 
150  do j = i1, i2
151  dothin = .false.
152  if( idx(j) .gt. 0) then
153  call kqsortr(x(1,j), si, idx(j))
154 ! call ksortinv(si, idx(j))
155  do i = 1, idx(j)-1
156  k1 = si(i)
157  if(k1 .lt. 0) goto 20
158  do l = i+1, idx(j)
159  k2 = si(l)
160  if(k2 .lt. 0) goto 10
161  if(abs(x(k1, j)- x(k2,j)) .gt. pixel ) goto 20
162 ! since xx is sorted, we may safe to skip
163 ! further check
164 !
165  if(abs(y(k1, j)- y(k2,j)) .gt. pixel ) goto 10
166  if(abs(z(k1, j)- z(k2,j)) .gt. pixel ) goto 10
167 ! remove k2
168  si(l) = - k2
169  dothin = .true.
170  10 continue
171  enddo
172  20 continue
173  enddo
174  if(dothin) then
175  nc0 = idx(j)
176  call mvdata(x, j, nc)
177  call mvdata(y, j, nc)
178  call mvdata(z, j, nc)
179  call mvdatai(codex, j, nc)
180  call mvdatai(chgx, j, nc)
181  idx(j) = nc
182 !/////////////
183  write(0,*) ' thinning ', nc0, "-->", nc
184 !//////////////
185  endif
186  endif
187  enddo
nodes z
nodes i
subroutine mvdata(f, j, nc)
Definition: timeslice.f:191
block data cblkElemag data *AnihiE ! Eposi< 1 TeV, anihilation considered *X0/365.667/, ! radiation length of air in kg/m2 *Ecrit/81.e-3/, ! critical energy of air in GeV *MaxComptonE/1./, ! compton is considered below 1 GeV *MaxPhotoE/1.e-3/, ! above this, PhotoElectric effect neg. *MinPhotoProdE/153.e-3/, ! below 153 MeV, no gp --> hadrons ! scattering const not MeV *Knockon true
Definition: cblkElemag.h:7
averg real MaxCPU integer idx(Maxp)
********************block data cblkHeavy ********************integer j data *HeavyG2symbol *data *HeavyG2code kiron data j
Definition: cblkHeavy.h:36
latitude latitude this system is used *****************************************************************! type coord sequence union map real y
Definition: Zcoord.h:25
block data cblkElemag data *AnihiE ! Eposi< 1 TeV, anihilation considered *X0/365.667/, ! radiation length of air in kg/m2 *Ecrit/81.e-3/, ! critical energy of air in GeV *MaxComptonE/1./, ! compton is considered below 1 GeV *MaxPhotoE/1.e-3/, ! above this, PhotoElectric effect neg. *MinPhotoProdE/153.e-3/, ! below 153 MeV, no gp --> hadrons ! scattering const not MeV *Knockon ! knockon is considered Obsolete *PhotoProd false
Definition: cblkElemag.h:7
real cut integer nc
Definition: Zprivate.h:1
subroutine mvdatai(f, j, nc)
Definition: timeslice.f:217
subroutine kqsortr(A, ORD, N)
Definition: kqsortr.f:23
! structure defining a particle at production ! Basic idea of what is to be contained in ! the particle structue is that dynamical ones should be included those derivable from the particle code ! is not included ******************************************************type fmom momentum sequence union map real e endmap map real * x
Definition: Zptcl.h:21
Here is the call graph for this function:
Here is the caller graph for this function:

◆ wtaskel()

subroutine wtaskel ( integer  fno,
real*4, dimension(np)  xa,
real*4, dimension(np)  ya,
real*4, dimension(np)  za,
integer*2, dimension(np)  code,
integer*2, dimension(np)  chg,
integer  np,
real*4  mulalpha 
)

Definition at line 359 of file timeslice.f.

References code2rgb().

Referenced by mkaskelfile(), and mkskelfiles().

359 !
360  implicit none
361  integer fno, np
362  real*4 xa(np), ya(np), za(np)
363  real*4 mulalpha
364  integer*2 code(np), chg(np)
365 
366  real r, g, b, alpha
367  integer i
368  if(np .gt. 0) then
369  write(fno, '(2i9)') np, np
370  do i = 1, np
371  write(fno, '(3g16.8)' ) xa(i), ya(i), za(i)
372  enddo
373  do i= 1, np
374  call code2rgb(code(i), chg(i), r, g, b, alpha)
375  write(fno, '("1 ", i8, 3f5.2,f6.3)')
376  * i-1, r, g, b, alpha*mulalpha
377  enddo
378  endif
*Zfirst p fm *Zfirst p Zfirst p Zfirst p *Zfirst p *Zfirst pos xyz r
Definition: ZavoidUnionMap.h:1
nodes i
*Zfirst p fm *Zfirst p Zfirst p code
Definition: ZavoidUnionMap.h:1
subroutine code2rgb(codei, chgi, r, g, b, alpha)
Definition: timeslice.f:381
dE dx *! Nuc Int sampling table g
Definition: cblkMuInt.h:130
real(4), save b
Definition: cNRLAtmos.f:21
integer, parameter fno
Definition: chook.f:43
Here is the call graph for this function:
Here is the caller graph for this function: