;----------------------------------------------------------------------------
;+
; NAME:
;       FORECAST
;
; PURPOSE:
;       The FORECAST procedure provides an interacting interface to:
;         - generate a simulated scan on the sky 
;         - display the scan over a dust map of the region
;         - display the scan over a CMB dipole map of the region
;         - calculate the dust contamination to the given scan
;           and display the result
;         - calculate the CMB dipole signal for the given scan 
;           and display the results 
;         - generate a Postscript file for hardcopy printing 
;
;       Presently, two types of scan are available:
;         - sinusoidal scan: the elevation of the beam 
;           is fixed, while the azimuth is modulated sinusoidally. 
;         - 360 degs scan: the elevation of the beam is fixed,
;           while the azimuth performs a full rotation on the sky. 
;
; CALLING SEQUENCE:
;       FORECAST, ra, dec, dust, dipole
;
; INPUTS/OUTPUTS:
;       ra: Right Ascension (hours)
;       dec: declination (degrees)
;       
;       These variables may have been generated in advance, or may be
;       generated interactively by the FORECAST procedure. 
;       In the former case they are considered as inputs,
;       otherwise they will contain in output the ra-dec generated by FORECAST.
;       
;       dipole: a vector containing the dipole signal in K for the 
;               given scan
;       dust: a vector containing the values of the 
;             equivalent CMB temperature fluctuation 
;             per sample due to dust, in K, for the given scan 
;
; PROCEDURE USED:
;       AZEL2RADEC: convert from horizon to celestial coordinates
;       DT_DUST: calculates dust contamination
;       DIPOLE: simulates dipole signal 
;       DIPOLE_GNOM: plot scan over dipole map - gnomic projection
;       DIPOLE_LAMBERT: plot scan over dipole map - lambert projection
;       DUST_GNOM: plot scan over dust map - gnomic projection
;       DUST_LAMBERT: plot scan over dust map - lambert projection
;       MAKE_SCAN: generates scan coordinates
;
; MODIFICATION HISTORY:
;       Created, Amedeo Balbi, August 1998
;       Modified, Amedeo Balbi, November 1998, to accept pre-generated
;       coordinates and to display the map in Lambert projection
;
;-
;----------------------------------------------------------------------------
pro FORECAST, ra, dec, dust, dipole

print,' '
print,' '
print,' ************** FORECAST V1.0 ************** '
print,' '
print,' '

print, '   (0)-generate scan'
print, '   (1)-use pre-generated scan'
read,igen

if igen eq 0 then MAKE_SCAN,time,az,el,ra,dec

print, '   (0)-gnomic projection (square patch smaller than 100X100 sq.degs)'
print, '   (1)-Lambert projection (full emisphere)'
read,iproj

; gnomic projection:
if iproj eq 0 then begin
    print, '   enter map side dimension (degrees) '
    print, '    ** maximum value: 100 degrees ** '
    read,side
center:
    print, '   (0)-automatic image center'
    print, '   (1)-manual image center'
    read,icenter
    if icenter then begin
        print, ' enter ra-dec (degrees, degrees) of image center:'
        read, racenter, deccenter
    endif else begin
        racenter=ra(n_elements(ra)/2)*15.
        deccenter=dec(n_elements(dec)/2)
    endelse

    loadct,5
    window,/free,xs=700,ys=700
    DUST_GNOM,racenter,deccenter,side,ra,dec
    wshow,iconic=0
    dust_image=tvrd()
    window,/free,xs=700,ys=700
    DIPOLE_GNOM,racenter,deccenter,side,ra,dec
    wshow,iconic=0
    dipole_image=tvrd()
    print, ' do you want to change the image center? '
    print, '   (0)-no'
    print, '   (1)-yes'
    read,ianswer
    if ianswer then goto, center
endif

; lambert projection
if iproj eq 1 then begin
    print, '   (0)-celestial coordinates '
    print, '   (1)-galactic coodinates'
    read,igal
    print, ' input +1 for northern emisphere, -1 for southern'
    read,ngsp
    loadct,5
    window,/free,xs=700,ys=700
    DUST_LAMBERT, ra, dec, ngsp, GAL_COORD=igal
    wshow,iconic=0
    dust_image=tvrd()
    window,/free,xs=700,ys=700
    DIPOLE_LAMBERT, ra, dec, ngsp, GAL_COORD=igal
    wshow,iconic=0
    dipole_image=tvrd()
endif

; working out dust contribution:
print, ' enter the frequency for dust extrapolation (GHz)'
print, '   (e.g. 150)'
read,nu
print, ' enter the spectral index for dust extrapolation'
print, '   (e.g. 1.6 (typical))'
read,alpha
dtitle=string('dust signal, frequency=',nu,' GHz, spectral index=',alpha)
nu=nu*1.e9                      ; converting to Hz
DT_DUST,ra,dec,nu,alpha,dust
window,/free
plot,dust,xtit='sample',ytit='dT_dust (K)',title=dtitle,psy=3

; working out dipole signal:
DIPOLE,ra,dec,dipole
window,/free
plot,dipole,xtit='sample',ytit='T_dipole (K)',title='dipole signal',psy=3

print,' do you want to generate a ps file with the results? '
print, '   (0)-no'
print, '   (1)-yes'
read,ips
if ips then begin
    psname=''
    print,' enter file name '
    read,psname
    psname1=psname+'1'
    psname2=psname+'2'
    psname3=psname+'3'
    print,''
    print,'   '+psname1+' will contain the dust map'
    print,'   '+psname2+' will contain the dipole map'
    print,'   '+psname3+' will contain the plot of the dust and dipole signals'
    set_plot,'ps'
    device,/portrait,xs=17.5,ys=16.7,yoff=5,file=psname1,/color,bits=8 
    tvscl,dust_image
    device,/close
    device,/portrait,xs=17.5,ys=16.7,yoff=5,file=psname2,/color,bits=8 
    tvscl,dipole_image
    device,/close
    !p.multi=[0,1,2]
    device,file=psname3,/portrait,xs=17.5,ys=21,yoff=3
    loadct,0
    mean=total(dust)/n_elements(dust)
    qdev=(dust-mean)^2
    var=total(qdev)/(n_elements(qdev)-1.)
    sdev=sqrt(var)
    stit=string('mean=',mean,' K, sdev= ',sdev,' K')
    plot,dust,xtit='sample',$
      ytit='dT_dust',title=dtitle,subtitle=stit,psy=3
    plot,dipole,xtit='sample',$
      ytit='T_dipole (K)',title='dipole signal',psy=3
    !p.multi=0
    device,/close
    set_plot,'x'
endif 

return
end