pro readall_res,dir=dir,res=res

;
; Name: READALL
;
; Called by: ANALYSIS.PRO
;  Or can be called on its own
;
; Program to read in all the McIntosh Archive (McA) Level 3 fits.gz files
; from directory dir, sample at one degree resolution latitude/longitude,
; and put into one big save file for later analysis.
;
; Written by Sarah Gibson -- updated September 2017
; modified by Tom Kuchar to include resultion keyword -- April 2020
;  Updated Nov 2022 to define feature numbering correctly in comments

; KEYWORD: DIR --  where LEVEL3FITS are
;          RES --  resolution to preserve in the archive file 
;                  the value determines if the resolution in the
;                  output latitude,longitude arrays are one degree
;                  (default) or higher, e.g. setting res=2 doubles the
;                  resolution of the array to 0.5 degreee
;

default,dir,'ARCHIVE/LEVEL3FITS/'

default,res,1

;
; to get the latest set of McA data, do the following:
;
; ftp -i ftp://ftp.ngdc.noaa.gov/STP/space-weather/solar-data/solar-imagery/composites/synoptic-maps/mc-intosh/ptmc_level3/ptmc_level3_fits/
; mget *
; quit
;
; Fits data include an image array that defines SOLAR MAGNETIC FEATURES
; for a particular Carrington map (naming convention includes CROT number).
; Also fits extensions contain arrays defining latitude/longitude.
; and also polarity of features (this allows plotting sunspot and plage separated by polarity)
; Each pixel in image array is assigned a number representing 
; a color which in turn defines the magnetic feature.
;


; Individual fits files can be used with programs "patmapcolortable.pro"
; to establish color table as below, and "plotfinal.pro" to make gifs
; (at higher resolution than one degree).
; These programs are provided at: 
;	https://www2.hao.ucar.edu/mcintosh-archive/plotting-software
;

;
; SOLAR MAGNETIC FEATURES
; each is identified by a single number in fits files
;

; (these next two color values should not exist in FITS files)
; 0 white
; 1 black

; (** indicates value in level3 FITS files
;
; **2 pos CH bndry
; **3 pos CH
; **4 pos polarity
; **5 neg CH bndry
; **6 neg CH
; **7 neg polarity
; **8 NL or PIL line
; **9 filaments
; **10 Sunspot (any polarity)
; **11 missing data
; **12 border
; **13 plage centers of changing flux outside active regions
; 14 neg sunspots (the FITS file won't have this, but reassigned below via polarity)
; 15 neg plage (the FITS file won't have this, but reassigned below via polarity)
;
; NOTE: The FITS file for later (e.g. WHPI) maps will have:
; **16 304 CH pos (284/211 later)
; **17 193 CH pos (195 EIT)
; **18 304 CH neg (284/211 later)
; **19 193 CH neg (195 EIT)
; **20 minimum aggregate coronal hole pos
; **21 minimum aggregate coronal hole neg
;

;
; Inputs: *fits.gz Level 3 McA files
;	sitting in directory dir
;
; Outputs: Save file named "mcintosharchive.dat" with:
;
; NLAT, NLON: 181, 361-- one degree resolution dimensions of Carrington maps
; NCROT: number of Carrington maps read in
; DATARRAY: the image array containing feature information
;	fltarr(ncrot,nlat,nlon) - varies in all 3 dimensions
; LATARRAY: latitude (-90=south pole, 90=north pole, 0=equator)
;	fltarr(ncrot,nlat,nlon) - varies in nlat dimension
; LONARRAY: longitude (0 - 360) 
;	fltarr(ncrot,nlat,nlon) - varies in nlon dimension
; POLARRAY: polarity (positive or negative)
;	fltarr(ncrot,nlat,nlon) - varies in all three dimensions
; CROTARRAY: Fractional Carrington longitude of central meridian passage
;		of each longitude of each map.
;		This is represented as a non-integer CROT number
;		where the non-integer value represents fraction
;		of full 360 degree longitude rotation. HOWEVER, because
;		sun rotates counterclockwise, longitude in the maps
;		decreases with increasing time, so that e.g. 
;		CROT XXXX.25=270 degrees, XXXX.75 = 90 degrees.
;	fltarr(ncrot,nlat,nlon) - varies in ncrot, nlon dimensions
; TIMEARRAY: Time corresponding to Carrington longitude at central meridian
;	fltarr(ncrot,nlat,nlon) - varies in ncrot, nlon dimensions
; CROTSMALLARRAY: Carrington longitude (latitude dimension collapsed)
;	fltarr(nlon*ncrot)
; TIMESMALLARRAY: time corresponding to CROTSMALLARRAY central meridian passage
;	fltarr(nlon*ncrot)
; CROTSMALLESTARRAY: Carrington rotation number for each map
;	fltarr(ncrot)
; TIMESMALLESTARRAY: time corresponding to beginning of each Carrington rotation
;	fltarr(ncrot)
;

;
; grab all the data, put into 3d cube
;  dimensions latitude, longitude, carrington rotation
;

 cd,dir,current=old_dir
 spawn,'ls -1 *fits.gz > archivelist'
 openr,1,'archivelist'
 fname=''
 allfiles=''
 count=0
 while ~ eof(1) do begin
  readf,1,fname
  if count gt 0 then allfiles=[allfiles,fname] else allfiles=fname
  count=count+1
 endwhile
 close,1
 spawn,'rm -f archivelist'

 ncrot=n_elements(allfiles)
 nlat=(180*res)+1
 nlon=(360*res)+1

 datarray=fltarr(ncrot,nlat,nlon)
 latarray=fltarr(ncrot,nlat,nlon)
 lonarray=fltarr(ncrot,nlat,nlon)
 polarray=fltarr(ncrot,nlat,nlon)
 crotarray=fltarr(ncrot,nlat,nlon)
 timearray=fltarr(ncrot,nlat,nlon)
 timesmallarray=fltarr(nlon*ncrot)
 timesmallestarray=fltarr(ncrot)

 index=0
 for i = 0,ncrot-1 do begin
     print,i,'  ', file_basename(allfiles(i))

;lonarray is in deg 0 to 360
     for j = 0,nlat-1 do begin
         lonarray(i,j,*)=findgen(nlon)/float(res)
     endfor
 
;latarray is in deg -90 to 90
     for k = 0,nlon-1 do begin
         latarray(i,*,k)=-90.+findgen(nlat)/float(res)
     endfor
 
; strip carrington rotation from filename
;   nameconvent='ptmc_compo_sm_'+yyyy+mm+dd+'_'+hh+mmin+ss+'_cr'+CR

  crot=strmid(allfiles(i),32,4)

  mreadfits,allfiles(i),header,image,/silent
  fits_read,allfiles(i),fullimage,exten_no=1
  fits_read,allfiles(i),longitude,exten_no=2
  fits_read,allfiles(i),latitude,exten_no=3
  fits_read,allfiles(i),polarity,exten_no=4

; higher resolution longitude/latutude indexing arrays
  rlongitude = fltarr((n_elements(longitude) -1)*res + 1)
  rlatitude  = fltarr((n_elements(latitude)  -1)*res + 1)

;  dlon = (longitude(n_elements(longitude) - 1) - longitude(0))/(n_elements(longitude) - 1)
;  dlat = (latitude(n_elements(latitude) - 1)   - latitude(0))/(n_elements(latitude) - 1)

  dlon = (shift(longitude,-1)-longitude)
  dlon(n_elements(longitude) - 1) = dlon(0)
  dlat = (shift(latitude,-1)-latitude)
  dlat(n_elements(latitude) - 1) = dlat(0)

  lonind= indgen(n_elements(longitude))
  latind= indgen(n_elements(latitude))

;  rlongitude((indgen(n_elements(longitude))*res))   = longitude
;  rlatitude((indgen(n_elements(latitude))*res))   = latitude

  for m=0,res-1 do begin
      rlongitude((lonind*res) + m) = longitude + (dlon/res)*m
      rlatitude((latind*res) + m)  = latitude + (dlat/res)*m
  endfor

; ensure the last elements are exact
  rlongitude(n_elements(rlongitude) - 1) = longitude(n_elements(longitude) - 1)
  rlatitude(n_elements(rlatitude) - 1) = latitude(n_elements(latitude) - 1)


; these remain the same regardless of resolution, lonmax not currently used
  if crot lt 1487 then begin
      lonmin=0
      lonmax=360
  endif else begin
      lonmin=30
      lonmax=390
  endelse

; modified here
; latitide and longitude are array elements of the larger image
; corresponding to one degree increments. Longitude array dimensions
; can be either 361 elements or 421 for later images, corresponding to
; the array element representing 0 to 360 deg or -30 to 390, respectively

; lonuse is always 0 to 360 deg, no matter what the dimensions of the
; original image

  if (res eq 1) then begin
      for j = 0,nlat-1 do begin
          for k = 0,nlon-1 do begin
              latuse=latitude(j)
              lonuse=longitude(k+lonmin)
              datarray(i,j,k)=fullimage(lonuse,latuse)
              polarray(i,j,k)=polarity(lonuse,latuse)
          endfor
      endfor
  endif else begin
      for j = 0,nlat-1 do begin
          for k = 0,nlon-1 do begin
              latuse=rlatitude(j)
              lonuse=rlongitude(k+lonmin*res)
              datarray(i,j,k)=fullimage(lonuse,latuse)
              polarray(i,j,k)=polarity(lonuse,latuse)
          endfor
      endfor
  endelse


;
; calculate time for each longitude, carrington rotation
;
  for k = 0,nlon-1 do begin
   timesmallarray(index)=anytim(carr2ex(crot+1.001-float(k)/360.))
   for j = 0,nlat-1 do timearray(i,j,k)=timesmallarray(index)
   index=index+1
  endfor

  if crot ne 1503 then begin
      timesmallestarray(i)=anytim(carr2ex(crot))
; this array is a float not a string
;  endif else timesmallestarray(i)='1966-01-09T09:36:08.64'
  endif else timesmallestarray(i)=-4.095014516420D+08

endfor

 crotarray(*,*,*)=(TIM2CARR(timearray(*,*,*),/dc))   
 crotsmallarray=(TIM2CARR(timesmallarray,/dc))   
 crotsmallestarray=(TIM2CARR(timesmallestarray,/dc))   

 cd,old_dir
 save,filename='mcintosharchive_res.dat',datarray,latarray,lonarray,polarray,crotarray,timearray,timesmallarray,timesmallestarray,nlon,nlat,ncrot,crotarray,crotsmallarray,crotsmallestarray
 
end