; PoCaFLIC.PRO
;+
; NAME:
;        POCAFLIC
; PURPOSE:
;       Creation of a file containing a series of Polar Camera
;       images transformed to global spherical coordinates at a
;       fixed height.  The user may select either PACE geomagnetic
;       coordinates or geographic coordinates.  The file is suitable
;        for display as a movie using SHOWFLIC.PRO.
; CATEGORY:
;
; CALLING SEQUENCE:
;        .RUN POCAFLIC
; INPUTS:
;        None.  The user is prompted for all needed inputs.
; KEYWORD PARAMETERS:
;        None.  The user is prompted for all needed inputs.
; OUTPUTS:
;        None.
; COMMON BLOCKS:
;       HEADBYTES: Contains all header bytes from the images read
;               in from the .RAY file, along with exposure start
;               UTs, exposure lengths, and flags indicating bad
;                images.  Shared with RDHBYTES.PRO.
;       COORDS: Contains various arrays and parameters used in
;               transforming the images to PACE coordinates, and
;               displaying the images.  Shared with MAP_PGM.PRO or
;                MAP_GEO.PRO.
; SIDE EFFECTS:
;       Creates a .PGM or .GEO file containing a sequence of 383x383
;       byte images, each one representing a single Polar Camera
;       image, mapped to PACE or geographic coordinates, respectively.
;       If the user wishes, the images are overlaid with a PACE or
;       geographic coordinate grid, and labelled with date, time,
;       and altitude information.  In the case of geographic
;       coordinates, coastlines are also plotted along with the
;        coordinate grid.
; RESTRICTIONS:
;
; PROCEDURE:
;
; MODIFICATION HISTORY:
;        Written 94/04 by D P Steele.  Documented 94/06 by D P Steele.
;        Modified 950316 by D P Steele: geographic mapping added as an
;                option.
;-

COMMON headbytes,oldfname,hbytes,utsec,exptimes,fwtemps,bad_image
COMMON coords,coordht,plat,plong,glat,glong,asize,az,za,col,row $
             ,vRlayer $
             ,eur_glat,eur_glong,eur_ht,eur_col,eur_row,eur_mlt $
             ,mlt_label_offset,plat_contours,plong_contours

t0=SYSTIME(1)

; load various string-valued arrays needed

@isitdos
take2=0B

; Determine graphics device to use.

dev=''
READ,'Enter desired graphics device ([Z]/X): ',dev
dev=STRUPCASE(dev)
IF dev NE 'X' THEN dev='Z'
!ORDER=0
IF !D.WINDOW EQ 0 THEN WDELETE,0
old_dev=!D.NAME
SET_PLOT,dev

; tell user what imaging modes and filters are available

cam_filt,filters

; prompt user and get input

cam=0 & filt=cam & day=cam & month=cam
PRINT,'Enter the year, month, day, camera, and filter desired'
READ,year,month,day,cam,filt

PRINT,'Enter the hour, minute, and second to start (default 00:00:00 UT)'
get_hms,sthms,shour,sminute,ssecond
IF sthms NE '' THEN start=[shour,sminute,ssecond] ELSE start=INTARR(3)

PRINT,'Enter the hour, minute, and second to finish (default 23:59:59 UT)'
get_hms,fihms,fhour,fminute,fsecond
IF fihms NE '' THEN finish=[fhour,fminute,fsecond] ELSE finish=[23,59,59]

eheight=[87,100,300,120,100,87,100,120,97,100]
PRINT,'Default emission height is '+STRTRIM(eheight(5*cam+filt),2) $
     +' km; enter preferred height or <RETURN> to accept this one:'
htstr=''
READ,htstr
IF htstr NE '' THEN eht=ROUND(FLOAT(htstr)) ELSE eht=eheight(5*cam+filt)

mtype=''
READ,'MLT/MLat dial or geographic coordinates? ([M]/G) ',mtype
mtype=STRUPCASE(mtype)
IF mtype EQ 'G' THEN BEGIN
    ctr=''
    READ,'Center the plot on the geomagnetic pole or Eureka? ([P]/E) ' $
       ,ctr
    pole=(STRUPCASE(ctr) NE 'E')
    IF NOT pole THEN BEGIN
        rotstr=''
        READ,'Rotate the plot to the correct "MLT" orientation? ([Y]/N) ' $
            ,rotstr
        rot=(STRUPCASE(rotstr) NE 'N')
        ENDIF ELSE rot=1B
    ENDIF ELSE BEGIN
    rot=0B
    pole=0B
    ENDELSE

resp=''
READ,'Do you want to label the image? ([Y]/N) ',resp
annotate=(STRUPCASE(resp) NE 'N')
CASE mtype OF
'M': BEGIN
        READ,'Do you want to show MLAT/MLT coordinates? ([Y]/N) ',resp
        grid=(STRUPCASE(resp) NE 'N')
        END
ELSE: BEGIN
        READ,"Show coastlines ('c') and/or coordinates ('g')? ",resp
        grid=STRUPCASE(resp)
        PRINT,'If you want to overplot a satellite footprint, enter'
        READ,'the standard 2-line element satellite name: ',resp
        sat=STRUPCASE(resp)
        END
    ENDCASE

IF (eht LT 300) AND (mtype NE 'G') THEN BEGIN
    READ,'Do you want to zoom in on the coordinate grid? ([Y]/N)',resp
    zoom=(STRUPCASE(resp) NE 'N')
    ENDIF ELSE zoom=0

READ,'Enter the van Rhijn background (R) to subtract: ',vRbg

; Prompt user for fixed saturation brightness

PRINT,'If you want to assign the top of the color bar to a fixed number of counts,'
PRINT,'enter that number now, otherwise enter anything else for auto-scaling.'
tops=''
READ,tops
flag=isnumber(tops,topnum)        ; topnum=0 if tops doesn't represent a number

; ensure a 4-digit year for display; make reasonable assumptions

IF year LT 1900 THEN BEGIN
    IF year LT 50 THEN year=year+2000
    IF year GT 50 THEN year=year+1900
ENDIF

; create title for image plot

site_names=['','Spyhill Airglow Observatory','Rabbit Lake, SK' $
           ,'Eureka AStrO Lab']
PCsite=PoCasite(year,month,day)
site=site_names(PCsite)

; Display only those pixels viewing the sky

getsig,year,month,day,mask0,mask1
CASE cam OF
 0: mask=mask0
 1: mask=mask1
    ENDCASE

; define path to desired images

template=filocate(year,month,day,cam,filt)

; Open the file as a series of small (header-size) units and get
; the header information

nimages=rdhbytes(template)

; Load (or create) the data to correct for the reduced transmission.

reltrans,oldfname,filtemp,filtrans

; Get the exposure time to be displayed

getexpt,exptimes,expt,selexp

; Select images to process

todo=getlims(nimages,utsec,start=start,finish=finish)

; Now load the colour table

cf=5*cam+filt
CASE cf OF
    2: ct=3
    7: ct=1
    8: ct=8
    ELSE: ct=0
    ENDCASE
loadct,ct,/silent

; Open the single big file of processed images and associate it to
; a variable.

OPENR,imunit,template,/GET_LUN
files=ASSOC(imunit,BYTARR(66048))

; Open another file to hold the transformed and annotated images.

IF mtype EQ 'G' THEN suffix='.geo' ELSE suffix='.pgm'
trname=pgmroot+dd $
      +STRING(year MOD 100,month,day,cam,filt,suffix $
             ,FORMAT='(3I2.2,2I1,A4)')
OPENW,trunit,trname,/GET_LUN
trans=ASSOC(trunit,BYTARR(383,383))

; Write information about the mapped image file contents to an
; appropriately-named information file.

IF mtype EQ 'G' $
    THEN geo_info,template,trname,cam,filt,year,month,day $
                 ,shour,sminute,ssecond,fhour,fminute,fsecond $
                 ,eht,vRbg,site,expt,topnum $
                 ,annotate=annotate,grid=grid,zoom=zoom,sat=sat $
    ELSE pgm_info,template,trname,cam,filt,year,month,day $
                 ,shour,sminute,ssecond,fhour,fminute,fsecond $
                 ,eht,vRbg,site,expt,topnum $
                 ,annotate=annotate,zoom=zoom

; Show how long it took to get this far

t1=SYSTIME(1)
PRINT,t1-t0,FORMAT='(F5.1," s to set up")'

; read in each image, decode the header, rescale from 8 bits,
; transform the image to the desired coordinates, and save it in
; the mapped image file.

ip=-1
FOR i=todo(0),todo(1) DO BEGIN
    whbad=WHERE(bad_image EQ i,nwhbad)
    IF nwhbad EQ 0 THEN BEGIN
        KIh=gethd(hbytes(*,i))

; check exposure time if necessary

        IF selexp THEN BEGIN
            thisexp=KIh.misc.exp_scale*KIh.exp.exposure/1000.
            go_on=(ABS(thisexp-expt) LT 1.0)
            ENDIF ELSE go_on=1

        IF go_on THEN BEGIN
            file=files(i)
            im=BYTE(file,512,256,256)
            im=cleanmod(im,3,2,3,/quiet)        ; fix hot pixels

; If file was renormalized to fit 8 bits, reverse this.

            IF file(290) GT 0 THEN BEGIN
                IF file(291) EQ 1 THEN im=LONG(im)*(2^file(290))
                IF file(291) EQ 0 THEN im=LONG(im)*file(290)
                ENDIF
; extract filter temperature, determine reduced filter transmission,
; and divide image by reduced transmission to recover actual emission
; brightnesses

            fwtemp=fwtemps(i)
            imtrans=interpol(filtrans,filtemp,fwtemp)
            im=ROUND(im/imtrans(0))

; Update the image pointer

            ip=ip+1

; Set the maximum displayed intensity

            IF topnum EQ 0 $
                THEN lmaxi=pcentile(im*mask,0.999) $
                ELSE lmaxi=topnum
            im(255,98)=lmaxi        ; force image to scale properly
            lmaxbh=pcentile(im*mask,0.001)

; Transform the image to PACE coordinates and label if desired

            IF mtype EQ 'G' $
                THEN map_geo,hbytes(*,i),im*mask,eht,pimg $
                            ,vRbg=vRbg,fullscale=[lmaxbh,lmaxi] $
                            ,zoom=zoom,annotate=annotate,grid=grid $
                            ,rot=rot,pole=pole,sat=sat $
                ELSE map_pgm,hbytes(*,i),im*mask,eht,88.664 $
                            ,-83.444,pimg $
                            ,vRbg=vRbg,fullscale=[lmaxbh,lmaxi] $
                            ,zoom=zoom,annotate=annotate

; Save the transformed image to the associated file

            trans(ip)=pimg
            EMPTY

; Update the information file with information about the current image

            pgm_imgs,trname,ip,hbytes(*,i),eur_mlt,plat(asize(1)/2,*)

; Show progress and check for user keypresses

            PRINT,'.',FORMAT='(A,$)'
            IF (ip+1) MOD 80 EQ 0 THEN PRINT,FORMAT='()'
            k=STRUPCASE(GET_KBRD(0))
            CASE 1 OF
      k EQ 'T': BEGIN                ; show progress and time of last image
                    PRINT,FORMAT='()'
                    PRINT,todo(0),i,todo(1),FORMAT='(3I4)'
                    PRINT,KIh.misc.tm
                    END
      k EQ 'X': GOTO,DONE        ; force exit RIGHT NOW
          ELSE:
                ENDCASE
            ENDIF
        ENDIF
    ENDFOR

DONE:
CLOSE,imunit & FREE_LUN,imunit
CLOSE,trunit & FREE_LUN,trunit
PRINT,FORMAT='()'
SET_PLOT,old_dev

t2=SYSTIME(1)
PRINT,t2-t1,ip+1,FORMAT='(F6.1," s to process ",I0," images")'

END