pro aurora_diff

  mydir='/Users/wzihan/rainbow'
  datestr='20130518'
  mapdir='/pngmap/'

  staname='sask'

  aele={i:0,j:0,elev:0.}
  elearr=replicate(aele,256,256)
  grd={lat:0.,lon:0.}
  grdarr=replicate(grd,257,257)
  cord={cenlat:0.,cenlon:0.,minlat:0.,minlon:0.,maxlat:0.,maxlon:0.}

  l=make_array(20,1,/integer,value=0)

  openr, 1, mydir+'/'+staname+'asiskymap.gm5' ; Skymap file made by me, containing the pixel AACGM coordinates using 2015 IRGF epoch
  readu, 1, elearr
  readu, 1, cord
  readu, 1, grdarr
  close, 1

  WINDOW,/FREE,XSIZE=480,YSIZE=480, retain=2
  xsize=!D.x_size
  ysize=!D.y_size

  image24 = BytArr(3, xsize, ysize)
  device, decomposed=0
  TVCRS, 0

  minlat=56.
  maxlat=64.5
  minlon=-55.
  maxlon=-35

  clat=(minlat+maxlat)/2.
  clon=(minlon+maxlon)/2.
  maxelv=10. ;
  restore, mydir+'/saskrainbow01.sav' ; rainbow image file 0430-0530 UT

  for hr=4,4 do begin
    for mn=30, 59 do begin
      for m=0,9 do begin

        loadct,0

        MAP_SET, /stere, clat,clon,LIMIT=[minlat,minlon,maxlat,maxlon],/noborder,charsize=2.

        for j=0,255 do begin
          for i=0,255 do begin
            if elearr[i,j].elev gt maxelv then begin
              ii=elearr[i,j].i
              jj=elearr[i,j].j
              mlat0=grdarr[jj,ii].lat
              mlat1=grdarr[jj,ii+1].lat
              mlat2=grdarr[jj+1,ii+1].lat
              mlat3=grdarr[jj+1,ii].lat
              mlon0=grdarr[jj,ii].lon
              mlon1=grdarr[jj,ii+1].lon
              mlon2=grdarr[jj+1,ii+1].lon
              mlon3=grdarr[jj+1,ii].lon

              if (m eq 0) and (mn eq 30) then begin
                trr1=0
                trg1=0
                trb1=0       
              endif else begin
                trr1=fix(rgb_struct_combined[(hr-4)*600+(mn-30)*10+m-1].red[ii,jj]*600.) >0
                trg1=fix(rgb_struct_combined[(hr-4)*600+(mn-30)*10+m-1].green[ii,jj]*600.) >0
                trb1=fix(rgb_struct_combined[(hr-4)*600+(mn-30)*10+m-1].blue[ii,jj]*600.) >0
              endelse
              
              trr2=fix(rgb_struct_combined[(hr-4)*600+(mn-30)*10+m].red[ii,jj]*600.) >0
              trg2=fix(rgb_struct_combined[(hr-4)*600+(mn-30)*10+m].green[ii,jj]*600.) >0
              trb2=fix(rgb_struct_combined[(hr-4)*600+(mn-30)*10+m].blue[ii,jj]*600.) >0

              ;chp2=trb2*65536L+trg2*256L+trr2
              diff=trg2-trg1
              if diff gt 127 then diff=127
              if diff lt -127 then diff=-127
              polyfill, [mlon0,mlon1,mlon2,mlon3],[mlat0,mlat1,mlat2,mlat3], color=diff+128
            endif
          endfor
        endfor

        map_grid, LATDEL=4.,LONDEL=10.,GLINETHICK=2.5,color='777777'x, latlab=minlon+.5,lonlab=minlat+.5,charsize=1.5,lats=60.,lons=0.,charthick=3.,/label
        tmstr=string(hr,mn,6*m,'(I02,I02,I02)')+' UT'
        xyouts, clon-4.,maxlat+.2,tmstr,color='FFFFFF'x,charsize=3.,charthick=4.

        Plots,-38.32,58.49,psym=7,symsize=1.3,color=250; C12
        xyouts, -38.32,58.49+0.1,'C12',color=250,charsize=1,charthick=1
        Plots,-51.24,60.64,psym=7,symsize=1.3,color=250; C06
        xyouts,-51.24,60.64+0.1,'C06',color=250,charsize=1,charthick=1

        image24=tvrd(/true)
        write_png, mydir+mapdir+tmstr+'_diff.png',image24
        ;cgColorbar, Divisions=4, Minor=5, Format='(F0.2)', Range=[-800, 800]
        endfor
    endfor
 endfor

  ;set_plot, 'win'
  print, 'finish'

end