pro readall,dir=dir

;
; 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
;

; KEYWORD: DIR --  where LEVEL3FITS are

default,dir,'ARCHIVE/LEVEL3FITS/'

;
; 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=181
 nlon=361

 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))

  for j = 0,nlat-1 do begin
   lonarray(i,j,*)=findgen(nlon)
  endfor

  for k = 0,nlon-1 do begin
   latarray(i,*,k)=-90.+findgen(nlat)
  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

   if crot lt 1487 then begin
    lonmin=0
    lonmax=360
   endif else begin
    lonmin=30
    lonmax=390
   endelse

  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
;
; 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.dat',datarray,latarray,lonarray,polarray,crotarray,timearray,timesmallarray,timesmallestarray,nlon,nlat,ncrot,crotarray,crotsmallarray,crotsmallestarray
 
end