Commit df9061f2 authored by Simone Iovenitti's avatar Simone Iovenitti
Browse files

Source code in IDL (need the astronomical library).

parent 2fe85cf6

star_coverage.pro

0 → 100644
+224 −0
Original line number Diff line number Diff line 
;+

; *****************************************

;              STAR COVERAGE

;

;

; Calcola la copertura angolare e il numero di stelle

; per un certo puntamento effettuato con il telescopio astri.

;       pcrab = queryvizier('USNO-B1',"m1",180, constraint='B1mag<8')

;       qcrab = queryvizier('USNO-B1',[83.629585, 22.01444444], 180, constraint='B1mag<8')

;         

;       [NB. SONO EQUIVALENTI QUESTE DUE!!! La seconda usa la posizione RA e DEC in GRADI di M1]

;

;

;-





    ROTATION = +1 ; +1 gira NATURAL, come il cielo (CCW), -1 gira in "REVERSE" (CW)

    VIEW     = +1 ; +1 mostra la SKYVIEW, -1 inserisce un flip N-S (CAMVIEW)   NB. Solo nei GRAFICI



;-- INSERISCI I DATI DELLA SIMULAZIONE:

    POINT = (catalogue("BRIGHT_STARS"))[where((catalogue("BRIGHT_STARS")).name eq "NCP")]  ;Navi ;Vega

    R     = 6.0   ; 5.8 ;deg

    MAG   = "5.5" ;"6.5"

    DATA  = '2019-02-26T02:36:43'   ; IN UTC!   ; DATA = sxpar( headfits(file) ,'DATE-OBS')

    LAST  = ten(2,30,00)            ; IN HMS!



;;-- PUNTAMENTO DI UN RUN DAL CATALOGO

;RUNS  = catalogue('RUNID') 

;RUN   = runs[0]

;POINT = {ra:run.ra,  dec:run.dec  }

;R     = 4.0 

;MAG   = "5.5" ;"6.5"

;DATA  = run.time_utc               ; IN UTC!   ; DATA = sxpar( headfits(file) ,'DATE-OBS')

;LAST  = ten(sec2hms(run.time))     ; IN HMS!

;tel   = telescope_parameters('ASTRI')                                                  ; Torna una struct con tutti i parametri del telescopio

;;********************



;;-- SORGENTI MINIARRAY

;a = indgen(15)  & remove,  [1, 2, 3, 4, 9, 10, 13], a  & print, a

;RUNS  = catalogue('MINIARRAY')

;R     = 4.5

;MAG   = "5.5" ;"6.5"

;DATA  = '2023-06-20T22:30:00'             ; IN UTC!   ; DATA = sxpar( headfits(file) ,'DATE-OBS')

;LAST  = ten(2,30,00)                      ; IN HMS!

;tel   = (telescope_parameters('MINIARRAY'))[1]                                                  ; Torna una struct con tutti i parametri del telescopio

;RUN   = runs[13]

;POINT = {ra:run.ra,  dec:run.dec  }

;;********************    

    

;;-- CONTROLLO CHE IL CAMPO SIA OSSERVABILE

    JD    = DATE_CONV(DATA, 'JD')

    EQ2HOR, POINT.RA, POINT.DEC, jd,                           alt_i, az_i, lat=tel.LAT_DEG,  lon=tel.lon_deg, altitude=tel.altitude

    EQ2HOR, POINT.RA, POINT.DEC, (jd*24.d + LAST*(0.5d))/24.d, alt_m, az_m, lat=tel.LAT_DEG,  lon=tel.lon_deg, altitude=tel.altitude

    EQ2HOR, POINT.RA, POINT.DEC, (jd*24.d + LAST)/24.d,        alt_f, az_f, lat=tel.LAT_DEG,  lon=tel.lon_deg, altitude=tel.altitude

    IF (ALT_I LT 10.0d) OR (ALT_M LT 10.0d) OR (ALT_F LT 10.0d) THEN BEGIN

      PRINT

      PRINT, "ATTENZIONE, IL CAMPO È BASSO SULL'ORIZZONTE!"

      PRINT

      ;STOP

    ENDIF

    

    

;;-- CONTROLLO CHE SIA SEMPRE NOTTE    

    sunpos, jd,                           ra_sun_i, dec_sun_i

    sunpos, (jd*24.d + LAST*(0.5d))/24.d, ra_sun_m, dec_sun_m

    sunpos, (jd*24.d + LAST)/24.d,        ra_sun_f, dec_sun_f

    EQ2HOR, ra_sun_i, dec_sun_i, jd,                           alt_sun_i, az_sun_i, lat=tel.LAT_DEG,  lon=tel.lon_deg, altitude=tel.altitude

    EQ2HOR, ra_sun_m, dec_sun_m, (jd*24.d + LAST*(0.5d))/24.d, alt_sun_m, az_sun_m, lat=tel.LAT_DEG,  lon=tel.lon_deg, altitude=tel.altitude

    EQ2HOR, ra_sun_f, dec_sun_f, (jd*24.d + LAST)/24.d,        alt_sun_f, az_sun_f, lat=tel.LAT_DEG,  lon=tel.lon_deg, altitude=tel.altitude

    IF (ALT_SUN_I GT -12.0d) OR (ALT_SUN_M GT -12.0d) OR (ALT_SUN_F GT -12.0d) THEN BEGIN

      PRINT

      PRINT, "ATTENZIONE, IL SOLE È QUASI ALL'ORIZZONTE."

      PRINT

      ;STOP

    ENDIF

    

    

    

;-- CONTROLLO LE POSIZIONI DI LUNA E PIANETI:

    mphase, jd, k   &   IF k GT 0.5 THEN PRINT, "     >> OCCHIO ALLA LUNA!"

    moonpos, jd, ram, decm, dism, /RADIAN 

    planet_coords, jd, rap, decp, /jd     ;'MERCURY','VENUS', 'MARS', 'JUPITER','SATURN','URANUS','NEPTUNE', 'MOON' ---> NB che la luna è sbagliata -.-

    ra_m_p = [ram, rap[0:-2]]

    dec_m_p = [decm, decp[0:-2]]

    gcirc,2, POINT.ra,POINT.dec, rap, decp, far   &  far = far/60./60.

    IF MIN(FAR) LE R THEN PRINT, "     >> OGGETTI DEL SISTEMA SOLARE NEL CAMPO DI VISTA!"

    

    

    

;-- ANGOLO PARALLATTICO (INIZIALE E SUA EVOLUZIONE):

    EQ2HOR, POINT.RA, POINT.DEC, jd, alt, az, ha, lat=tel.LAT_DEG,  lon=tel.lon_deg, altitude=tel.altitude   ; Trovo HA ; ALTAZ2HADEC, alt, az, tel.lat_deg, ha, dec_on_date  ; HA = LST - RA  ma fai la precess! (LST:   CT2LST, lst, tel.lon_deg, 1,ten(20,36), 26, 02, 2019 )

    parang_i  = PARANGLE(ha, POINT.DEC, tel.lat_deg, /DEGREE)                             ; HA cresce linearmente col tempo, PARANG invece no!

    parang_m  = PARANGLE(ha+(LAST*(1./2.)*15.), POINT.DEC, tel.lat_deg, /DEGREE)          ; QUesto intermedio serve per riconoscere il verso di rotazione...

    parang_f  = PARANGLE(ha+(LAST*15.), POINT.DEC, tel.lat_deg, /DEGREE)                  ; 

    parang    = -angoli_ordinati(parang_i, parang_m, parang_f)                            ; *** METTO IL MENO PER FARLO GIRARE CCW COME IN GENERE FANNO GLI ANGOLI ***





    

;-- CERCO LE STELLE NEL FOV, NELLA BANDA SPETTRALE GIUSTA:

    JUMP=0 & if nel(lista) ne 0 then begin & jump=1 & goto, here & endif                           ; vuol dire che stai runnando più volte, forse un debug?

    lista = queryvizier('USNO-B1',[POINT.RA, POINT.DEC], R*60.d, constraint='B1mag<'+MAG) ; radius in arcmin, B1=Optical B band between 400 and 500 nm

    IF nel(lista) eq 0 THEN BEGIN

       IF lista[0] eq -1 THEN BEGIN & print

            print, "LISTA VUOTA!"

            perimetro = bordo_camera("REDUCED")

            plot_stars_i = plot([0.], [0.], linestyle=6, symbol=3, sym_filled=0, sym_size=8, sym_color='black', ASPECT_RATIO=1, XTITLE='[deg]', YTITLE='[deg]',  TITLE="ASTRI FOV", xrange=[perimetro.x/tel.PlateScaledeg-1, perimetro.x/tel.PlateScaledeg+1], yrange=[perimetro.x/tel.PlateScaledeg-1, perimetro.x/tel.PlateScaledeg+1])

            plot_bordo   = plot(perimetro.x/tel.PlateScaledeg, perimetro.y/tel.PlateScaledeg, color='light blue', /OVERPLOT)

            e = ELLIPSE(0., 0., MAJOR=dms2deg(0,22,0), /DATA, FILL_COLOR="black", THICK=0., TRANSPARENCY=50., PATTERN_ORIENTATION=45, PATTERN_SPACING=7, PATTERN_THICK=3)

            print & stop

       ENDIF

    ENDIF

    

    field = []

    

    FOR ii=0, nel(lista)-1 DO BEGIN    

      querysimbad, "USNO-B1.0 "+lista[ii].USNO_B1_0,  ra_simbad, dec_simbad, id_simbad, found=found

      

      IF FOUND THEN nome = id_simbad

      IF NOT FOUND THEN BEGIN                                                             ; per alcune stelle SIMBAD non consulta USNO [-.-...], e io guardo le stelle che ci sono su Hipparcos!

        hip = queryvizier("I/239/hip_main", [lista[ii].raj2000,lista[ii].dej2000], 3)     ; arcmin. Se c'è da fare la PRECESSIONE usa: precess, ra, de, 2000, 1991.25, /PRINT      

        distanze = []

        for jj=0, nel(hip)-1 do distanze = [distanze, deg(acos( sin(rad(lista[ii].DEJ2000))*sin(rad(hip[jj]._DE_ICRS)) + cos(rad(lista[ii].DEJ2000))*cos(rad(hip[jj]._DE_ICRS))*cos(rad( lista[ii].RAJ2000 -hip[jj]._RA_ICRS)) )) ]

        nome = "HIP" +  string(hip[where(distanze eq min(distanze))].hip, FORMAT='(I0)') 

      ENDIF

      

      far = deg(acos( sin(rad(lista[ii].DEJ2000))*sin(rad(POINT.dec)) + cos(rad(lista[ii].DEJ2000))*cos(rad(POINT.dec))*cos(rad( lista[ii].RAJ2000 -POINT.ra)) ))

      ;EQUIVALENTE: gcirc,2, POINT.ra,POINT.dec, lista[3].RAJ2000, lista[3].DEJ2000, far   &  far = far/60./60.

      posang, 1, deg2hrs(POINT.ra), POINT.dec, deg2hrs(lista[ii].RAJ2000), lista[ii].DEJ2000, ang      ; LA SORG 1 È IL CENTRO! decimal hours - decimal degrees - deg    

      cirrange, ang            ; *** NB. Metto i posang da 0 a 360 (quando plotti l'angolo zero ce l'hai a destra) Poi sommo (90-parang_iniziale) per avere il FoV orientato come ASTRI. ***   

      

      field = [field, {name:nome,     ra:lista[ii].RAJ2000,    de:lista[ii].DEJ2000,     mag:lista[ii].B1MAG,      distance:far,    posang:double(ang + 90.d - parang_i)    } ]

    ENDFOR

    

    

    

;-- STAMPO IL CAMPO STELLARE INQUADRATO CON I NOMI:

    r0 = double( ceil(max(field.distance)+1))                           ; Per verifica su Stellrium:

    angoli = [signum(VIEW)*rad(field.posang)]                           ; deg2hms(field[4].ra, /string)

    starplot = plot(field.distance*cos(angoli), field.distance*sin(angoli), linestyle=6, symbol=3, sym_filled=1, sym_size=8, sym_color='red', ASPECT_RATIO=1, XTITLE='[deg]', YTITLE='[deg]', xrange=[-r0-1, r0+1], yrange=[-r0-1, r0+1], TITLE="NOMI STELLE INQUADRATE")

    t = text(field.distance*cos(angoli), field.distance*sin(angoli), string(indgen(nel(field)), FORMAT='(I0)') + "-" +field.name,  /DATA)





    here:

    if JUMP EQ 1 THEN BEGIN & print & print, "LA LISTA ERA PRE-COMPILATA" & print & ENDIF





;-- STAMPO I RISULTATI

    step   = 0.01d ; 36''

    center = dms2deg(0,22,0)

    

    copertura = angular_coverage(field, parang=SIGNUM(ROTATION)*parang, step=step, center=center)

    perimetro = bordo_camera("REDUCED")

    

    print, nel(copertura.covered)

    print, nel(where(copertura.covered eq 1.d))

    

    

    ;-- calcolo le x e le y per fare i plot (e per usare fuori_camera che lavora in x,y)

    passi     = FINDGEN((abs(parang)/step)+1, INCREMENT = SIGNUM(ROTATION)*signum(parang)*step) 

    stars_i_x = field.distance*cos(angoli)

    stars_i_y = field.distance*sin(angoli)

    stars_f_x = field.distance*cos(angoli+signum(VIEW)*SIGNUM(ROTATION)*rad(parang))

    stars_f_y = field.distance*sin(angoli+signum(VIEW)*SIGNUM(ROTATION)*rad(parang))

    trails_x  = [] & for kk=0, nel(field)-1 do trails_x = [trails_x, field[kk].distance*cos( signum(VIEW)*rad(field[kk].posang+passi) ) ]

    trails_y  = [] & for kk=0, nel(field)-1 do trails_y = [trails_y, field[kk].distance*sin( signum(VIEW)*rad(field[kk].posang+passi) ) ]

    coverd_x  = r0*cos( signum(VIEW)*rad(copertura.angles[where(copertura.covered eq 1)]))

    coverd_y  = r0*sin( signum(VIEW)*rad(copertura.angles[where(copertura.covered eq 1)]))

    

    TITOLO = "SKY-VIEW" & IF (VIEW EQ -1) THEN TITOLO = "CAM-VIEW"

    plot_covered = plot(coverd_x,  coverd_y,  linestyle=6, symbol=3, sym_filled=1, sym_size=1, sym_color='red', xrange=[-r0-1, r0+1], yrange=[-r0-1, r0+1], ASPECT_RATIO=1, TITLE="ASTRI FOV IN "+TITOLO)

    plot_stars_i = plot(stars_i_x, stars_i_y, linestyle=6, symbol=3, sym_filled=0, sym_size=8, sym_color='black', ASPECT_RATIO=1, XTITLE='[deg]', YTITLE='[deg]', xrange=[-r0-1, r0+1], yrange=[-r0-1, r0+1], /overplot)

    plot_stars_f = plot(stars_f_x, stars_f_y, linestyle=6, symbol=3, sym_filled=1, sym_size=8, sym_color='black', ASPECT_RATIO=1, XTITLE='[deg]', YTITLE='[deg]', xrange=[-r0-1, r0+1], yrange=[-r0-1, r0+1], /overplot)

    plot_trails  = plot(trails_x,  trails_y,  linestyle=6, symbol=3, sym_filled=1, sym_size=1, sym_color='black', /overplot)

    plot_bordo   = plot(perimetro.x/tel.PlateScaledeg, perimetro.y/tel.PlateScaledeg, color='light blue', /OVERPLOT)

    e = ELLIPSE(0., 0., MAJOR=center, /DATA, FILL_COLOR="black", THICK=0., TRANSPARENCY=50., PATTERN_ORIENTATION=45, PATTERN_SPACING=7, PATTERN_THICK=3)

    t = text(0.05, 0.8, "COPERTURA:")

    t = text(0.05, 0.75,string(double(nel(coverd_x))/(360./step)*100., FORMAT= '(D0.2)') + "%."   )

    t = text(0.05, 0.7, "FUORI CAMERA.")

    print

    print, "Fuori dalla camera la copertura è al ", string(double(nel(coverd_x))/(360./step)*100., FORMAT= '(D0.2)') ,"%."  

    print









;-- Controllo i FUORI_CAMERA

    

    ro = double(ceil(max(sqrt(perimetro.x^2.d + perimetro.y^2.d))/tel.PlateScaledeg))

    

    FUORI_CAMERA, stars_i_x*tel.PlateScaledeg, stars_i_y*tel.PlateScaledeg, stars_i_x_c, stars_i_y_c

    FUORI_CAMERA, stars_f_x*tel.PlateScaledeg, stars_f_y*tel.PlateScaledeg, stars_f_x_c, stars_f_y_c

    FUORI_CAMERA, trails_x*tel.PlateScaledeg,  trails_y*tel.PlateScaledeg,  trails_x_c,  trails_y_c

    

    trails_x_c_ok = trails_x_c[where(sqrt(trails_x_c^2.d + trails_y_c^2.d)/tel.PlateScaledeg gt center)]

    trails_y_c_ok = trails_y_c[where(sqrt(trails_x_c^2.d + trails_y_c^2.d)/tel.PlateScaledeg gt center)]

    

    plot_stars_i = plot(stars_i_x_c/tel.PlateScaledeg, stars_i_y_c/tel.PlateScaledeg, linestyle=6, symbol=3, sym_filled=0, sym_size=8, sym_color='black', ASPECT_RATIO=1, XTITLE='[deg]', YTITLE='[deg]',  TITLE="ASTRI FOV IN "+TITOLO, xrange=[-ro-1, ro+1], yrange=[-ro-1, ro+1])

    plot_stars_f = plot(stars_f_x_c/tel.PlateScaledeg, stars_f_y_c/tel.PlateScaledeg, linestyle=6, symbol=3, sym_filled=1, sym_size=8, sym_color='black', /overplot)

    plot_trails  = plot(trails_x_c_ok/tel.PlateScaledeg,   trails_y_c_ok/tel.PlateScaledeg, linestyle=6, symbol=3, sym_filled=1, sym_size=1, sym_color='black', /overplot)

    plot_trails  = plot(trails_x_c/tel.PlateScaledeg,   trails_y_c/tel.PlateScaledeg, linestyle=6, symbol=3, sym_filled=1, sym_size=1, sym_color='black', /overplot)

    plot_bordo   = plot(perimetro.x/tel.PlateScaledeg, perimetro.y/tel.PlateScaledeg, color='light blue', /OVERPLOT)

    e = ELLIPSE(0., 0., MAJOR=center, /DATA, FILL_COLOR="black", THICK=0., TRANSPARENCY=50., PATTERN_ORIENTATION=45, PATTERN_SPACING=7, PATTERN_THICK=3)



;-- Disegno il tratto rosso per il coverage

    test = round(cirrange(deg(atan(trails_y_c_ok, trails_x_c_ok)))*100.d)/100.d    ; NB. Passare da test ad angola è necessario non per il grafico (andrebbe bene anche con test) 

    angola = DINDGEN((360.d/step) , INCREMENT=step, START=0.)                ;     ma per fare il calcolo corretto delle percentuali!

    match, angola, test, suba, subb, COUNT = count                ;angola[suba] è uguale test[subb] e contiene tutti i valori in comune tra i due array



    camera_cover_x = ro*cos(signum(VIEW)*rad(angola[suba]))

    camera_cover_y = ro*sin(signum(VIEW)*rad(angola[suba]))

    plot_covered = plot(camera_cover_x,  camera_cover_y,  linestyle=6, symbol=3, sym_filled=1, sym_size=1, sym_color='red', /OVERPLOT)

    t = text(0.05, 0.8, "COPERTURA:")

    t = text(0.05, 0.75,string(double(count)/nel(angola)*100., FORMAT= '(D0.2)') +"%." )    

    t = text(0.05, 0.7, "NELLA CAMERA.")

    

    

;-- Stampo i risultati a video:  

    print

    print, "Fuori dalla camera la copertura è al ", string(double(nel(coverd_x))/(360./step)*100., FORMAT= '(D0.2)') ,"%."

    print, "Dentro la camera la copertura è al ", string(double(count)/nel(angola)*100., FORMAT= '(D0.2)') ,"%."

    print, "L'evoluzione dell'angolo parallattico vale: ", string(parang, FORMAT='(D0.1)')

    print

    

END