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
|
!
! aperture photometry tester
!
! Copyright © 2019-20 F.Hroch (hroch@physics.muni.cz)
!
! This file is part of Munipack.
!
! Munipack is free software: you can redistribute it and/or modify
! it under the terms of the GNU General Public License as published by
! the Free Software Foundation, either version 3 of the License, or
! (at your option) any later version.
!
! Munipack is distributed in the hope that it will be useful,
! but WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
! GNU General Public License for more details.
!
! You should have received a copy of the GNU General Public License
! along with Munipack. If not, see <http://www.gnu.org/licenses/>.
!
program testaphot
use fitsaphot
use daofotometr
implicit none
integer :: naper, i
real :: r0, e, f, g, g2
real, dimension(:), allocatable :: raper, apsum
character(len=666) :: arg
call get_command_argument(1,arg)
read(arg,*) r0
call get_command_argument(2,arg)
read(arg,*) e
call rfits('art.fits',naper,raper,apsum)
do i = 1, naper
f = 1 - exp(-raper(i)**2/2/r0**2)
g = apsum(i) / apsum(naper)
g2 = g2d(raper(i),r0,e)
write(*,'(6f10.3)') raper(i),g,f,g-f,g2,g2-f
! grow(raper(i)/r0)*grow(raper(i)/r0/sqrt(1-e**2))
end do
deallocate(raper,apsum)
contains
real function grow(r)
real, intent(in) :: r
grow = (1 + erf(r/sqrt(2.0))) / 2
end function grow
real function g2d(raper,r0,e)
real, intent(in) :: raper,r0,e
integer, parameter :: ndim = 2500
! real, dimension(-ndim:ndim,-ndim,ndim) :: g
integer :: n,i,j
real :: s,c,g,r2
n = 0
s = 0
c = 100
do i = -ndim, ndim
do j = -ndim, ndim
r2 = (i**2 + j**2) / c**2
! if( .true. ) then
if( r2 <= raper**2 ) then
g = exp(-r2/r0**2/2)
s = s + g
n = n + 1
end if
end do
end do
! write(*,*) s,n,s / n / c,s/250591.047
g2d = s / n / c
g2d = s/250591.047
g2d = s / (c*ndim)
end function g2d
subroutine rfits(filename,naper,raper,apsum)
use titsio
character(len=*), intent(in) :: filename
integer, intent(out) :: naper
real, dimension(:), allocatable, intent(out) :: raper, apsum
character(len=80) :: key
real, dimension(:), allocatable :: col
integer :: status,nrows,irow,i
type(fitsfiles) :: fits
logical :: anyf
irow = 1
status = 0
call fits_open_table(fits,filename,FITS_READONLY,status)
call fits_get_num_rows(fits,nrows,status)
allocate(col(nrows))
call fits_read_key(fits,'NAPER',naper,status)
allocate(raper(naper),apsum(naper))
do i = 1, naper
write(key,'(a,i0)') 'APER',i
call fits_read_key(fits,key,raper(i),status)
end do
do i = 1, naper
call fits_read_col(fits,3+2*i,1,0.0,col,anyf,status)
apsum(i) = col(irow)
end do
deallocate(col)
call fits_close_file(fits,status)
call fits_report_error(error_unit,status)
end subroutine rfits
end program testaphot
|
