pro hexpoint
;program to calculate the centres for a hexagonal 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,outdat,"hexpoint.ann",/get_lun
openw,listlun,"hexpoint.cat",/get_lun
printf,outdat,"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))*sin(60./!RADEG)*bw/2.
if not nyquist then begin
	dec_row_incr=dec_row_incr*2.
	ra_incr=ra_incr*2.
endif

decposns0=dec+(findgen(nr)+1-(nr+1)/2)*dec_row_incr
decposns1=decposns0(0:nr-2)+dec_row_incr/2.
rapos=ra+(findgen(nc)+1)*ra_incr-ra_incr*(nc+1)/2.

for z=0,nc-1 do begin
	if z/2. mod 1 eq 0 then for k=0,nr-1 do $
	printf,outdat,"circle w ",rapos(z),decposns0(k),bw/2.
	if z/2. mod 1 ne 0 then for k=0,nr-2 do $
	printf,outdat,"circle w ",rapos(z),decposns1(k),bw/2.
endfor

rapos_str=strarr(nc)
decpos0_str=strarr(nr)
decpos1_str=strarr(nr-1)

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 i=0,nc-1 do begin
ra_in_dec=dec2deg(rapos(i),/ra)
ra_in_dec=[string(fix(ra_in_dec[0:1])),string(format='(f5.2)',ra_in_dec(2))]
for j=0,2 do if ra_in_dec(j) lt 10 then $
ra_in_dec(j)=strcompress(string(0)+ra_in_dec(j),/remove_all)
rapos_str(i)=strcompress(string(format='(a2,a2,a)',$
ra_in_dec(0),ra_in_dec(1),ra_in_dec(2)),/remove_all)
endfor

for i=0,nr-1 do begin
dec_in_dec=dec2deg(decposns0(i),/dec)
dec_in_dec=[string(fix(dec_in_dec[0:1])),string(format='(f5.2)',dec_in_dec(2))]
for j=0,2 do if dec_in_dec(j) lt 10 then $
dec_in_dec(j)=strcompress(string(0)+string(dec_in_dec(j)),/remove_all)
decpos0_str(i)=strcompress(dec_in_dec(0)+dec_in_dec(1)+dec_in_dec(2),/remove_all)
endfor

for i=0,nr-2 do begin
dec_in_dec=dec2deg(decposns1(i),/dec)
dec_in_dec=[string(fix(dec_in_dec[0:1])),string(format='(f5.2)',dec_in_dec(2))]
for j=0,2 do if dec_in_dec(j) lt 10 then $
dec_in_dec(j)=strcompress(string(0)+dec_in_dec(j),/remove_all)
decpos1_str(i)=strcompress(dec_in_dec(0)+dec_in_dec(1)+dec_in_dec(2),/remove_all)
endfor

for z=0,nc-1 do begin
 for k=0,nr-1 do begin
  name_str=strcompress("point"+string(z+1)+"_"+string(k+1),/remove_all)+" "
  if z/2. mod 1 eq 0 then begin 
   printf,listlun,name_str,rapos_str(z)," ",decpos0_str(k),endstr
  endif else begin
   if k ne nr-1 then printf,listlun,name_str,rapos_str(z)," ",decpos1_str(k),endstr 
  endelse
 endfor
endfor





npoints=(floor(nc/2)+1)*nr+ceil(nc/2)*(nr-1)

print,strcompress("Finished making a "+string(fix(nc))+'x'+string(fix(nr))+" hexagonal pointing plan")
if nyquist then print,strcompress("Sampled at NYQIST RATE with "+string(fix(npoints))+" points.")$
else print,strcompress("Sampled at BEAM width with "+string(fix(npoints))+" points.")

endprog:
close,/all
end