; history
;
;  2012-09-29 written
;  2012-10-01 output to blesklog.txt changed to mill.txt
;  2012-10-06 bug in insertion of data gaps fixed in mill.ps
;  2012-10-13  vertical scale redefined as E^(1/3)to show fair waether field,
;              min(blue) and max(red) lines added to the average (black)
;  2013-03-13 thick black line, thin dotted red and blue lines
;  2013-03-18 thick black line, thin solid red and blue lines
;  2013-06-16 5min detail on the right
;  2013-06-18 bugs fixed
;  2013-06-20 proportions of overview and detail changed
;
function ver
  return,'2013-06-20'
end

pro processmill
;
  print, '---- Electric field mill v. '+ver()
  f=file_search('*.csv',count=nf)
  if nf gt 0 then begin
    f=f[sort(f)]
    for i=0,nf-1 do begin
      readmill,f[i], ts,kvm
      akvm=total(kvm)/n_elements(kvm)
      minkvm=min(kvm,max=maxkvm)
      addlog,f[i],ts[0],ts[n_elements(ts)-1],akvm
      addfig,f[i],ts[0],ts[n_elements(ts)-1],akvm,minkvm,maxkvm, ts,kvm
    end
  end
  print, '---- Done.'
end


pro readmill,file, ts,kvm
;
  print, '----  Reading mill file '+file
     openr,1,file & a='' & r=0
     while not eof(1) do begin
       readf,1,a & r=r+1
       if r eq 1 then ts=strmid(a,0,26) else ts=[ts,strmid(a,0,26)]
       v=0. & s=0 & reads,strmid(a,27,strlen(a)-27),v,s
       if r eq 1 then kvm=v else kvm=[kvm,v]
     end
     close,1
  print, '        -> Reading OK '+strmid(ts[0],0,19)+'-'+strmid(ts[n_elements(ts)-1],11,8)
end

pro addlog,file, ts0,ts1,akvm
    openw,2,'mill.txt'
    printf,2,'   '+strmid(ts0,11,8)+'-'+strmid(ts1,11,8)+$
                 string(akvm,form='(f7.2)')
    close,2
end

pro addfig,file, ts0,ts1,akvm,minkvm,maxkvm, ts,kvm
;__accumulated overview figure-----------------
   maxtim=4ull*86400000000000ull ; 4 days,
   if keyword_set(ts0) then str_tns,ts0+'.000.000',tns0
   if keyword_set(ts1) then str_tns,ts1+'.000.000',tns1
   ;__ data need to exist
   prev= file_test('mill_data.idl',/read) and $
             not file_test('mill_data.idl',/zero)
   ;__ data need to be a vector
   if prev then begin ;-------
       restore,'mill_data.idl'
       siz=size(data)
       if siz[0] eq 1 then begin
         pw=0.33
         data=[[data*(1-0.03/data^pw)^(1./pw)],[data],[data*(1+0.03/data^pw)^(1./pw)]]
         siz=size(data)
       end
       if siz[0] ne 2 then prev=0
       if n_elements(tim) eq 0 then prev=0
       if siz[1] ne n_elements(tim) then prev=0
       if siz[2] ne 3 then prev=0
   end
    ;__ no backward steps in time
   if prev then begin
      if keyword_set(ts0) then tim=[tim,tns0]
      ntim=n_elements(tim)
      outt=0
      while tim[ntim-1]-tim[outt] gt maxtim do outt=outt+1
      tim=tim[outt:ntim-1] & ntim=ntim-outt
      if ntim lt 2 then prev=0
    end
   if prev then begin
      datan=fltarr(ntim,3)
      datan[0:ntim-2,*]=data[outt:ntim+outt-2,*]
      if keyword_set(ts0) then datan[ntim-1,*]=[[minkvm],[akvm],[maxkvm]] else datan[ntim-1,*]=data[ntim+outt-1,*]
      data=datan
      ;
     dtim=long64(tim[1:ntim-1]) - long64(tim[0:ntim-2])
     ibad=where(dtim le 0 ,nbad)
     if nbad gt 0 then prev=0
   end
   if prev then begin
      medtim=median(dtim,/doub)
      igap=where(dtim gt 2*medtim,ngap)
      if ngap gt 0 then begin
        gap=lonarr(ngap+1)
        for ig=1,ngap do gap[ig]=round(dtim[igap[ig-1]]/medtim)
        igap=[-1,igap,ntim-1]
      end else begin
        gap=[0]
        igap=[-1,ntim-1]
      end
      ;
      datan=replicate(!values.f_nan,siz[1]+1+total(gap),3)
      for ig=0,ngap do $
        datan[igap[ig]+1+total(gap[0:ig]):igap[ig+1]+total(gap[0:ig]),*]=$
         data[igap[ig]+1:igap[ig+1],*]
      if keyword_set(ts0) then print,'----    Added value '+string(akvm,form='(f6.2)')+' at '+strmid(ts0,0,19)+' to the mill data.'
   end else begin
       print,'----    First value '+string(akvm,form='(f6.2)')+' to the mill data at '+strmid(ts0,0,19)
       data=[[minkvm],[akvm],[maxkvm]] & tim=tns0
   end
   save,file='mill_data.idl',tim,data
   millfile=strmid(file,0,15)+'_mill.ps'
   plot_init,millfile,yfrac=0.3
   !p.charsize=0.65
   if not keyword_set(ts0) then begin
     tns1=tim[ntim-1]
     tns_str,tns1,ts0 & ts1=ts0
   end

   timeh=rotate(-findgen(n_elements(datan[*,0])),2)/float(n_elements(datan[*,0])-1)*$
            float((tns1[0]-tim[0])/3600000000000ull)
   pw=0.33
   ticks=[20,10,5,2,0.7,0.1]
   tickn=['20','10','5','2','0.7','0.1']
   ;
   ; overview
   iok=where(finite(datan),nok)
   if nok gt 0 then begin
     ine=where(datan[iok] lt 0.,nne)
     ipo=where(datan[iok] ge 0.,npo)
     if npo gt 0 then datan[iok[ipo]]=datan[iok[ipo]]^pw
     if nne gt 0 then datan[iok[ine]]=-abs(datan[iok[ine]])^pw
   end
   tns_str,tim[0],tss
   plot,timeh, datan[*,1],xtit='Time (hours) from '+strmid(ts0,0,19), $
         ytit='Vertical electric field (kV/m)',col=255,yticks=n_elements(ticks)*2,$
         tit='Overview plot from '+strmid(tss,0,19)+' to '+strmid(ts1,0,19),$
         ytickv=[-(ticks^pw),rotate(ticks^pw,2)],$
         ytickname=['-'+tickn,rotate(tickn,2)],$
         yticklen=0.01,xticklen=0.05,xstyle=1,yran=[-25^pw,25^pw],/nodata,pos=[0.045,0.14,0.64,0.92]
    oplot,[timeh[0],0],-[0,0],line=0,col=255
    for i=0,n_elements(ticks)-1 do begin
      oplot,[timeh[0],0],[ticks[i],ticks[i]]^pw,line=1,col=255
      oplot,[timeh[0],0],-[ticks[i],ticks[i]]^pw,line=1,col=255
    end
    oplot,timeh, datan[*,0],psym=-3,col=64,thick=1
    oplot,timeh, datan[*,2],psym=-3,col=1,thick=1
    oplot,timeh, datan[*,1],psym=-3,col=255,thick=5
    ;
    oplot,[timeh[0]*0.95,timeh[0]*0.9], [-17^pw,-17^pw],psym=-3,col=255,thick=5
    xyouts,timeh[0]*0.9,-19^pw,'Average', col=255,/data
    oplot,[timeh[0]*0.75,timeh[0]*0.7], [-17^pw,-17^pw],psym=-3,col=1,thick=1
    xyouts,timeh[0]*0.7,-19^pw,'Maximum', col=1,/data
    oplot,[timeh[0]*0.55,timeh[0]*0.5], [-17^pw,-17^pw],psym=-3,col=64,thick=1
    xyouts,timeh[0]*0.5,-19^pw,'Minimum', col=64,/data
    ;
    ; last 6.5 min
   iok=where(finite(kvm),nok)
   if nok gt 0 then begin
     ine=where(kvm[iok] lt 0.,nne)
     ipo=where(kvm[iok] ge 0.,npo)
     if npo gt 0 then kvm[iok[ipo]]=kvm[iok[ipo]]^pw
     if nne gt 0 then kvm[iok[ine]]=-abs(kvm[iok[ine]])^pw
   end
   tss=ts[0]
   nts=n_elements(ts)
   ts1=ts[nts-1]
   str_tns,ts+'.000.000',tns
   timem=(tns-tns[0])/60d9
   plot,timem, kvm,xtit='Time (minutes) from '+strmid(tss,0,19), $
         ytit='Vertical electric field (kV/m)',col=255,yticks=n_elements(ticks)*2,$
         tit='Last interval: '+strmid(tss,0,19)+' - '+strmid(ts1,11,8),$
         ytickv=[-(ticks^pw),rotate(ticks^pw,2)],$
         ytickname=['-'+tickn,rotate(tickn,2)],$
         yticklen=0.01,xticklen=0.05,xstyle=1,yran=[-25^pw,25^pw],/nodata, /noerase, pos=[0.70,0.14,0.999,0.92]
    oplot,[0,timem[nts-1]],-[0,0],line=0,col=255
    for i=0,n_elements(ticks)-1 do begin
      oplot,[0,timem[nts-1]],[ticks[i],ticks[i]]^pw,line=1,col=255
      oplot,[0,timem[nts-1]],-[ticks[i],ticks[i]]^pw,line=1,col=255
    end
    oplot,timem, kvm,psym=-3,col=255,thick=3
   ;
   plot_close
   file_copy,millfile,'mill.ps',/overwrite
end





;---------------------------------------------------------------
;--------     general time routines, nanosecond MJD    ---------
;---------------------------------------------------------------



function into_tns,nanosecond,second,minute,hour,day,month,year
 ;
 ; composition of modified jul. date in nanoseconds,
 ;   8byte unsigned integer, nanoseconds from 17nov1858 till 6jun2443
 ;
  return,nanosecond +1000000000ull* (second+60ull*(minute+60ull*(hour + $
                  24ull*(julday(Month, Day, Year)- 2400001ull) ) ) )
end


pro from_tns, tns, nanosecond,second,minute,hour,day,month,year
 ;
 ; decomposition of modified jul. date in nanoseconds,
 ;   8byte unsigned integer, nanoseconds from 17nov1858 till 6jun2443
 ;
  CALDAT, long(tns    /86400000000000ull +2400001ull), Month , Day, Year
  nanosecond=ulong( tns mod 1000000000ull)
  second    =byte( (tns/    1000000000ull) mod 60ull)
  minute    =byte( (tns/   60000000000ull) mod 60ull)
  hour      =byte( (tns/ 3600000000000ull) mod 24ull)
end


pro str_tns, str,tns
 ;
 ; string -> modified jul. date in nanoseconds in 8byte unsigned integer
 ;
  n=n_elements(str)
  Year=uintarr(n) & Month=bytarr(n) & Day=bytarr(n)
  Hour=bytarr(n) & Minute=bytarr(n) & Second=bytarr(n) & nanosecond=ulonarr(n)
  Yeari=0ul &  Monthi=0ul &  Dayi=0ul &  Houri=0ul &  Minutei=0ul
  Secondi=0ul &  millisecondi=0ul &  microsecondi=0ul & nanosecondi=0ul
  for i=0,n-1 do begin
     reads,str(i),Yeari, Monthi, Dayi, Houri, Minutei, Secondi, millisecondi, microsecondi,nanosecondi, $
          form='(i4,1x,i2,1x,i2,1x,i2,1x,i2,1x,i2,1x,i3,1x,i3,1x,i3)'
     Year(i)=Yeari & Month(i)=Monthi & Day(i)=Dayi
     Hour(i)=Houri & Minute(i)=Minutei & Second(i)= Secondi
     nanosecond(i)=nanosecondi+1000ul*microsecondi+1000000ul*millisecondi
  end
  tns=into_tns(nanosecond,second,minute,hour,day,month,year)
end


pro tns_str,  tns,str
 ;
 ;  modified jul. date in nanoseconds in 8byte unsigned integer -> string
 ;
   from_tns, tns, nanosecond,second,minute,hour,day,month,year
  n=n_elements(tns)
  str=strarr(n)
  for i=0l,n-1 do $
    str(i)=string(Year(i), Month(i), Day(i), Hour(i), Minute(i), Second(i), $
                  nanosecond(i)/1000000ul, (nanosecond(i) mod 1000000ul)/1000ul,$
                  nanosecond(i) mod 1000ul,$
         form='(i4.4,"-",i2.2,"-",i2.2,"T",i2.2,":",i2.2,":",i2.2,".",i3.3,".",i3.3,".",i3.3)')
end





;---------------------------------------------------------------
;--------------     general plotting routines    ---------------
;---------------------------------------------------------------




pro plot_init,file,yfrac=yfrac
 ;
 ; open lanscape ps plot
 ; close by device,/close
 ;
 set_plot,'ps'
 if not keyword_set(yfrac) then yfrac=1.
 device,file=file,xsize=25.,ysize=18.*yfrac,/landscape,xoffset=1.,yoffset=28.,$
             bits_per_pixel=8,/color, /helvetica
 !p.font=0
 !p.charsize=1.
 !p.thick=1
 !x.thick=2
 !y.thick=2
 bw=0
 color_palette, bw, rp,gp,bp, colors
 tvlct,rp,gp,bp
end


pro plot_close
;
; closes ps file, restores normal crt plotting
 !p.multi=0
 device,/close
 if !VERSION.OS_FAMILY eq 'Windows' then set_plot,'win' else set_plot,'x'
end





pro color_palette, bw, rp,gp,bp, colors
;+
; color_palette
;   defines  a color palette and several particular colors
; Nov, 29, 1999, O.Santolik, written
; Jun, 06, 2000, O.Santolik, bw scale corrected
;-
 if bw then begin
   rp=[255,round(findgen(63)/62.*250.),255-round(findgen(190)/189.*255.),85,0]
   gp=rp & bp=rp
 endif else begin
  rp1=[255,255,250,246,242,238,234,230,226,222,218,214,210,206,202,198,194,190,$
       186,182,178,174,170,165,161,157,153,149,145,141,137,133,129,125,121,117,$
       113,109,105,101, 97, 93, 89, 85, 80, 76, 72, 68, 64, 60, 56, 52, 48, 44,$
    40, 36, 32, 28, 24, 20, 16, 12,  8,  4,  0,  0,  0,  0,  0,  0,  0,  0,$
     0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0]
  rp2=[  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,$
     0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,$
     0,  0,  4,  8, 12, 16, 20, 24, 28, 32, 36, 40, 44, 48, 52, 56, 60, 64,$
    68, 72, 76, 80, 85, 89, 93, 97,101,105,109,113,117,121,125,129,133,137,$
       141,145,149,153,157,161,165,170,174,178,182,186,190,194,198,202,206,210]
  rp3=[214,218,222,226,230,234,238,242,246,250,255,255,255,255,255,255,255,255,$
       255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,$
       255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,$
       255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,$
       255,255,128,  0]
  rp=[rp1,rp2,rp3]
  gp1=[255,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,$
     0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,$
     0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,$
     0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  7, 15, 23, 31, 39, 47, 55,$
    63, 71, 79, 87, 95,103,111,119,127,135,143,151,159,167,175,183,191,199]
  gp2=[207,215,223,231,239,247,255,255,255,255,255,255,255,255,255,255,255,255,$
       255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,$
       255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,$
       255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,$
       255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255]
  gp3=[255,255,255,255,255,255,255,255,255,255,255,250,246,242,238,234,230,226,$
       222,218,214,210,206,202,198,194,190,186,182,178,174,170,165,161,157,153,$
       149,145,141,137,133,129,125,121,117,113,109,105,101, 97, 93, 89, 85, 80,$
    76, 72, 68, 64, 60, 56, 52, 48, 44, 40, 36, 32, 28, 24, 20, 16, 12,  8,$
     4,  0,128,  0]
  gp=[gp1,gp2,gp3]
  bp1=[255,  0,  4,  8, 12, 16, 20, 24, 28, 32, 36, 40, 44, 48, 52, 56, 60, 64,$
    68, 72, 76, 80, 85, 89, 93, 97,101,105,109,113,117,121,125,129,133,137,$
       141,145,149,153,157,161,165,170,174,178,182,186,190,194,198,202,206,210,$
       214,218,222,226,230,234,238,242,246,250,255,255,255,255,255,255,255,255,$
       255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255]
  bp2=[255,255,255,255,255,255,255,246,238,230,222,213,205,197,189,180,172,164,$
       156,148,139,131,123,115,106, 98, 90, 82, 74, 65, 57, 49, 41, 32, 24, 16,$
     8,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,$
     0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,$
     0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0]
  bp3=[  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,$
     0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,$
     0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,$
     0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,$
     0,  0,128,  0]
  bp=[bp1,bp2,bp3]
 endelse
 colors={white:0,black:255,grey:254,red:1,blue:64,$
     green:125,yellow:190,cyan:96,violet:33,orange:220}
end
