aboutsummaryrefslogtreecommitdiff
path: root/pkg/tbtables/fitsio/ftgsvd.f
blob: c7e1d30bb48269069178e09e3daa79b393d22089 (plain) (blame)
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
C----------------------------------------------------------------------------
        subroutine ftgsvd(iunit,colnum,naxis,naxes,blc,trc,inc,
     &  nulval,array,anynul,status)

C       read a subsection of double precision data values from an image or
C       a table column.  

C       iunit   i  fortran unit number
C       colnum  i  number of the column to read from
C       naxis   i  number of dimensions in the FITS array
C       naxes   i  size of each dimension. 
C       blc     i  'bottom left corner' of the subsection to be read
C       trc     i  'top right corner' of the subsection to be read
C       inc     i  increment to be applied in each dimension
C       nulval  i  value that undefined pixels will be set to
C       array   i  array of data values that are read from the FITS file
C       anynul  l  set to .true. if any of the returned values are undefined
C       status  i  output error status
C
C       written by Wm Pence, HEASARC/GSFC, June 1993

        integer iunit,colnum,naxis,naxes(*),blc(*),trc(*),inc(*),status
        double precision array(*),nulval
        logical anynul,anyf

C-------COMMON BLOCK DEFINITIONS:--------------------------------------------
        integer nb,ne
        parameter (nb = 20)
        parameter (ne = 200)
        integer bufnum,chdu,hdutyp,maxhdu,hdstrt,hdend,nxthdr,dtstrt
        integer nxtfld
        logical wrmode
        common/ft0001/bufnum(199),chdu(nb),hdutyp(nb),maxhdu(nb),
     &  wrmode(nb),hdstrt(nb,ne),hdend(nb),nxthdr(nb),dtstrt(nb),nxtfld
C-------END OF COMMON BLOCK DEFINITIONS:------- -----------------------------

        integer i,i1,i2,i3,i4,i5,i6,i7,i8,i9,row,rstr,rstp,rinc
        integer str(9),stp(9),incr(9),dsize(10)
        integer felem,nelem,nultyp,ninc,ibuff,numcol
        logical ldummy
        character caxis*20

C       this routine is set up to handle a maximum of nine dimensions

        if (status .gt. 0)return

        if (naxis .lt. 1 .or. naxis .gt. 9)then
                status=320
                write(caxis,1001)naxis
1001            format(i20)
                call ftpmsg('NAXIS ='//caxis//' in the call to FTGSVD '
     &          //'is illegal.')
                return
        end if

C       if this is a primary array, then the input COLNUM parameter should
C       be interpreted as the row number, and we will alway read the image
C       data from column 2 (any group parameters are in column 1).

        ibuff=bufnum(iunit)
        if (hdutyp(ibuff) .eq. 0)then
C               this is a primary array, or image extension
                if (colnum .eq. 0)then
                    rstr=1
                    rstp=1
                else
                    rstr=colnum
                    rstp=colnum
                end if
                rinc=1
                numcol=2
        else
C               this is a table, so the row info is in the (naxis+1) elements
                rstr=blc(naxis+1)
                rstp=trc(naxis+1)
                rinc=inc(naxis+1)
                numcol=colnum
        end if

        nultyp=1
        anynul=.false.
        i1=1
        do 5 i=1,9
                str(i)=1
                stp(i)=1
                incr(i)=1
                dsize(i)=1
5       continue
        do 10 i=1,naxis
                if (trc(i) .lt. blc(i))then
                        status=321
                        write(caxis,1001)i
        call ftpmsg('In FTGSVD, the range specified for axis '//
     &  caxis(19:20)//' has the start greater than the end.')
                        return
                end if
                str(i)=blc(i)
                stp(i)=trc(i)
                incr(i)=inc(i)
                dsize(i+1)=dsize(i)*naxes(i)
10      continue

        if (naxis .eq. 1 .and. naxes(1) .eq. 1)then
C               This is not a vector column, so read all the rows at once
                nelem=(rstp-rstr)/rinc+1
                ninc=rinc
                rstp=rstr
        else
C               have to read each row individually, in all dimensions
                nelem=(stp(1)-str(1))/inc(1)+1
                ninc=incr(1)
        end if

        do 100 row=rstr,rstp,rinc
         do 90 i9=str(9),stp(9),incr(9)
          do 80 i8=str(8),stp(8),incr(8)
           do 70 i7=str(7),stp(7),incr(7)
            do 60 i6=str(6),stp(6),incr(6)
             do 50 i5=str(5),stp(5),incr(5)
              do 40 i4=str(4),stp(4),incr(4)
               do 30 i3=str(3),stp(3),incr(3)
                do 20 i2=str(2),stp(2),incr(2)

        felem=str(1)+(i2-1)*dsize(2)+(i3-1)*dsize(3)+(i4-1)*dsize(4)
     &  +(i5-1)*dsize(5)+(i6-1)*dsize(6)+(i7-1)*dsize(7)
     &  +(i8-1)*dsize(8)+(i9-1)*dsize(9)

        call ftgcld(iunit,numcol,row,felem,nelem,ninc,
     &  nultyp,nulval,array(i1),ldummy,anyf,status)
        if (status .gt. 0)return
        if (anyf)anynul=.true.
        i1=i1+nelem

20              continue
30             continue
40            continue
50           continue
60          continue
70         continue
80        continue
90       continue
100     continue
        end