pro gridpoint
;program to calculate the centres for a square pointing pattern
;indicate the number the ra and dec of the first pointing, and the 
;number of pointings in the centre row, and the numbe of rows
;program also generates a text file.

;HISTORY
; Created feb24 Erik M.

read,"Coords of centre pointing: hh mm ss.s sdd mm ss.s: ",rah,ram,ras,decd,decm,decs
read,"Enter velocity (km/s): ",vel
vframe=strarr(1)
read,"Enter velocity type (i.e. helio, tropo, geo, LSR: h/t/g/l): ",vframe
if total(strlowcase(vframe(0)) eq ['h','t','g','l']) eq 0 then begin
    print,'Not understood'
    goto,endprog
endif

vtype=strarr(1)
read,"Enter velocity type (i.e. optical, z(radio), z(optical): o/zr/zo): ",vtype
if total(strlowcase(vtype(0)) eq ['o','zr','zo']) eq 0 then begin
    print,'Not understood'
    goto,endprog
endif

read,"Number of pointings in centre row, number of rows: ", nc,nr
read,"Enter beamwidth in minutes: ",bw
nrate=strarr(1)
read,"Sample at Nyquist rate? (y/n): ", nrate
if strlowcase(nrate) eq 'y' then nyquist=1 else nyquist=0
bw=bw/60.

openw,outann,"gridpoint.ann",/get_lun
openw,outcat,"gridpoint.cat",/get_lun
rastring=strarr(3)
decstring=strarr(3)
printf,outann,"colour red"

ra=deg2dec(rah,ram,ras,/ra)
dec=deg2dec(decd,decm,decs,/dec)
dec_row_incr=bw/2.
ra_incr=(1/cos(dec/!RADEG))*bw/2.
if not keyword_set(nyquist) then begin
	dec_row_incr=dec_row_incr*2.
	ra_incr=ra_incr*2.
endif
decposns0=(findgen(nr)-(nr)/2)*dec_row_incr+dec


ra=ra-(ra_incr*(nc-1)/2.)
rastr=strarr(3)
decstr=strarr(3)
case strlowcase(vframe(0)) of
'h':vframe='helio'
't':vframe='topo'
'g':vframe='geo'
'l':vframe='lsr'
endcase

case strlowcase(vtype(0)) of
'o': vtype='optical'
'zr': vtype='z(radio)'
'zo': vtype='z(optical)'
endcase
endstr=strcompress(" j "+string(vel)+" "+vframe+" "+vtype)

for z=0,nc-1 do begin
	for k=0,nr-1 do begin
            printf,outann,"circle w ",ra,decposns0(k),bw/2.
            ra2=ra
            dec2=decposns0(k)
            rastr=string(dec2deg(ra2,/ra))
            decstr=string(dec2deg(dec2,/dec))
            rastr[0:1]=string([fix(rastr(0)),fix(rastr(1))])
            decstr[0:1]=string([fix(decstr(0)),fix(decstr(1))])
            rastr(2)=string(format='(f5.2)',rastr(2))
            decstr(2)=string(format='(f5.2)',decstr(2))
   if rastr(0) lt 10 then rastr(0)=strcompress(string(0)+rastr(0),/remove_all)
   if rastr(1) lt 10 then rastr(1)=strcompress(string(0)+rastr(1),/remove_all)
   if rastr(2) lt 10 then rastr(2)=strcompress(string(0)+rastr(2),/remove_all)
if decstr(0) lt 10 then decstr(0)=strcompress(string(0)+decstr(0),/remove_all)
if decstr(1) lt 10 then decstr(1)=strcompress(string(0)+decstr(1),/remove_all)
if decstr(2) lt 10 then decstr(2)=strcompress(string(0)+decstr(2),/remove_all)
rastr=strjoin(strcompress(rastr,/remove_all),"")
decstr=strjoin(strcompress(decstr,/remove_all),"")
vel=string(vel)
titlestr=strcompress("point_"+string(z+1)+"_"+string(k+1),/remove_all)

printf,outcat,strcompress(titlestr+" "+rastr+" "+decstr+endstr)

        endfor
ra=ra+ra_incr
endfor

print,strcompress("Finished making a "+string(fix(nc))+'x'+string(fix(nr))+" square grid pointing plan.")

endprog:
close,/all
end