; ID. plotsun.pro idl 5.5 18/9/2002
;+
; NAME:
;   PlotSun
;
; PURPOSE:
;   widget to compute and plot solar path for a given day, year,
;   longitude and latitude
;
; CATEGORY:
;   Astronomy
;   solar energy
;
; CALLING SEQUENCE:
;
;   PlotSun
;
; INPUTS:
;   from the input window:
;   latitude: degrees and decimal degrees
;   longitude: degrees and decimal degrees
;   time zone
;   time interva in minutes for computation
;   year
;   day of the year (1 - 365)
;   output file name
;     Default values are for Chamonix
;
; OUTPUTS:
;
;   graphic output : polar diagram of sun's path
;   and ASCII file with azimuth and elevation angles
;
; SIDE EFFECTS:
;    Creates a new window on every run
;
; REQUIREMENTS
;   functions solarpos, julianday, polardiagram
;     help/info files plotsun.hlp, infoJD.hlp
;
; MODIFICATION HISTORY:
;   Javier Corripio,  jgc@geo.ed.ac.uk
;   Created:  August, 2001
;   Modified:     September, 2002
;
;-
;  Copyright (c) 2002 Javier G Corripio.
;  jgc@geo.ed.ac.uk  http://www.geo.ed.ac.uk/~jgc/
;  This software may be used, copied, or redistributed as long as it is not
;  sold and this copyright notice is reproduced on each copy made.  This
;  routine is provided as is without any express or implied warranties
;  whatsoever.

@solarpos
@julianday
@polardiagram


;	Function circle is from Paul Krummel:
; 	http://cow.physics.wisc.edu/~craigm/idl/archive/msg00655.html
; 		http://www.dar.csiro.au/
;;       From his routine PLOT_WIND_ROSE, which inspired some bits of this
;	  routine
;
; MODIFICATION HISTORY:
; 	Written by:	Javier G. Corripio September, 2002
;
;-

FUNCTION CIRCLE, xcenter, ycenter, radius, nPts=nPts

IF N_Elements(nPts) EQ 0 THEN nPts = 100
points = (2 * !PI / (nPts-1)) * FIndGen(nPts)
x = xcenter + radius * Cos(points)
y = ycenter + radius * Sin(points)
RETURN, Transpose([[x],[y]])

END

PRO POLARDIAGRAM,radius,theta,time,title

	;; rotate angles to position 0 at North
thetarot=!pi/2-theta*!dtor
	;; radius is given as elevations, change to zenithal angle
radius_z = 90.-radius

PLOT, radius_z, thetarot, /nodata, $
        /polar, /isotropic, $
        xstyle=5, ystyle=5, $
        yrange=[-90,90],xrange=[-90,90],$
      	XMARGIN=[15,5],YMARGIN=[4,6]
        ;TITLE='polar diagram of solar path',$

bg=replicate(1,640)#indgen(512)
tv,(BYTSCL(bg)-50 > 0)
loadct,1

XYOUTS, .5,.95,title,ALIGNMENT=.5,/norm,charsize=1.0

	;; plot zenithal concentric circles @10 degrees interval
zenit_circ=findgen(10)*10/90.
FOR i=0,10-1 do begin
        PlotS, circle(0, 0, zenit_circ[i]*90, nPts=361)
	xyouts,0, zenit_circ[i]*90+1,' '+$
		strtrim(string(round(zenit_circ[i]*90)),2),color='cacaca'x
	xyouts,0, -zenit_circ[i]*90+1,' '+$
		strtrim(string(round(zenit_circ[i]*90)),2),color='cacaca'x
ENDFOR

	;; Plot 'nr' radial lines at 360/nr degree interval
nr=24
anglenr=2*!pi/float(nr)
rdl_r = replicate(90,nr)  ;; 90 is the maximum zenithal angle
rdl_th = findgen(nr)*anglenr
FOR i=0,nr-1  DO OPLOT, [0,rdl_r[i]],[0,rdl_th[i]],/Polar

	;; plot the labels for the radial lines
l_theta=STRTRIM(STRING(ROUND(rdl_th*!RADEG)),2)

rdl_polar = FLTARR(2,nr)
rdl_polar[0,*] = !PI/2.-rdl_th	; makes 0 at the top
rdl_polar[1,*] = rdl_r + 5	; 5 points offset from the outer circle border

	;; label in rectangular coordinates (value of angle)
l_theta_rect = CV_COORD(from_polar=rdl_polar,/to_rect)

XYOUTS, l_theta_rect[0,*],l_theta_rect[1,*],l_theta,ALIGNMENT=0.5

;; Ytitle
XYOUTS, 0.05,.5,'zenith angle',ALIGNMENT=.5,/norm,orientation=90,$
	charsize=1.0,color='cacaca'x

	;; plot the values
	colours
;OPLOT, radius_z,thetarot, /polar, thick=3, color=52479L
OPLOT, radius_z,thetarot, /polar, thick=3, color=!clr.gold

	;; plot the times at hh:00
oclock=where(ABS(FIX(time)-time) LT 1e-3,count)
IF count EQ 0 THEN oclock = 0
OPLOT, radius_z[oclock],thetarot[oclock], /polar, PSYM=1, color=65535L
hours_label = strtrim(string(fix(time[oclock])),2)+'h'
hours_pos_polar=FLTARR(2,size(oclock,/dim))
hours_pos_polar[0,*]=thetarot[oclock]
hours_pos_polar[1,*]=radius_z[oclock]+10
hours_pos_rect = CV_COORD(from_polar=hours_pos_polar,/to_rect)
XYOUTS, hours_pos_rect[0,*],hours_pos_rect[1,*],hours_label,ALIGNMENT=.5,$
	color=52479L, charsize=1.0

END

PRO plotsun_event, ev

COMMON plotsun_block, latitude,longitude,timezone,yearw,julianday,interval,$
    sundata,outfname

  widget_control, ev.id, get_uvalue=value
  if (n_elements(value) eq 0) then value = ''
  name= widget_info(ev.id,/uname)

case (name) of
"DONE": WIDGET_CONTROL, /destroy, ev.top

"HELP":begin
    cd,CURRENT=pwd
    XDISPLAYFILE, FILEPATH('plotsun.hlp',ROOT_DIR=pwd), $
    GROUP = ev.top, TITLE = 'Solar Position Help'
    end

"COMPUTEJD": begin
        JulianDay
       end
"COMPUTE": begin
    WIDGET_CONTROL, latitude, GET_VALUE = latit
    WIDGET_CONTROL, longitude, GET_VALUE = longit
    WIDGET_CONTROL, timezone, GET_VALUE = timez
    WIDGET_CONTROL, yearw, GET_VALUE = yearval
    WIDGET_CONTROL, julianday, GET_VALUE = julD
    WIDGET_CONTROL, interval, GET_VALUE = intv
    WIDGET_CONTROL, outfname, GET_VALUE = outfnme

    lat = FLOAT(latit[0])
    longi = FLOAT(longit[0])
    tmz = FLOAT(timez[0])
    year = FIX(yearval[0])
    JD = FIX(julD[0])
    interv = intv[0]/60.0
    fname = STRING(outfnme[0])


sundata={year:year,JDay:JD,time:0.00,timezone:tmz,latitude:lat,longitude:longi}

wnum=!D.window

sp=solarpos(sundata)
sunrise=fix(sp.sunrise)
sunset=fix(sp.sunset)+1
sundata.time=sunrise
declination = sp.declination

polarday=0
IF (ABS(lat) GT 66.5) and (lat/ABS(lat) NE declination/ABS(declination)) THEN BEGIN
    window,wnum+1
    XYOUTS, .5,.5,align=.5,charsize=1.5,/norm,$
       'day and latitude correspond to POLAR NIGHT'
    RETURN
ENDIF ELSE BEGIN
;;polar day
    IF ABS(TAN(lat*!DTOR)*TAN(declination*!DTOR)) GT 1  THEN BEGIN
       sunrise=0.0
       sunset=24.0
       polarday=1
    ENDIF
ENDELSE

timesteps=ROUND( ((sunset-sunrise)/interv) )+2

sol=findgen(3,timesteps)*0.0

for i=0,timesteps-1 do begin
    sp=solarpos(sundata)
    sol(0,i)=sundata.time
    sol(1,i)=sp.elevation
    sol(2,i)=sp.azimuth
    sundata.time=sundata.time+interv
endfor

;*** write output file
OPENW,2,fname
latlongstring=strcompress(string(FIX(lat)))+'!Uo!N latitude,'+$
    strcompress(string(FIX(longi)))+'!Uo!N longitude'

JDlong=julday(1,0,year)+JD
CALDAT,JDLong,M,D,Y

printf,2,'Solar position for the '+strcompress(string(D))+' /'+$
    strcompress(string(M))+' /'+strcompress(string(Y))
printf,2,' at  '+strcompress(string(lat))+' latitude, '+ $
       strcompress(string(longi))+' longitude

printf,2,FORMAT='(A5,2A14)', 'time','elevation','azimuth'
printf,2, FORMAT='(A5,2A14)','...','...','...'

timeplot=FLTARR(timesteps)
FOR i=0,timesteps-1 do begin
    ;**convert decimal fraction of hours to minutes
    hour=FIX(sol[0,i])
    min=(sol[0,i]-hour)*60/100.0
    hourmin=hour+min
    IF (.60-min LT 0.01) THEN hourmin = hour + 1 ;;avoids hh:60min
    timeplot[i] = hourmin
    PRINTF,2,FORMAT='(F5.2,2F)',hourmin,sol[1,i],sol[2,i]
ENDFOR
CLOSE,2
    cd,CURRENT=pwd
    XDISPLAYFILE, FILEPATH(fname,ROOT_DIR=pwd), $
    GROUP = ev.top, TITLE = 'Computed solar elevation and azimuth'

;;plot results
;device, decomposed = 0
loadct, 0
wnum=!D.window
window,wnum+1,xsize=600
title = 'Solar path for the '+strcompress(string(D))+' /'+$
    strcompress(string(M))+' /'+strcompress(string(Y))+$
    ' at  '+latlongstring

polardiagram, sol[1,*],sol[2,*],timeplot,title

;; compute sunrise and sunset of the upper limb of the sun
;;  ;; this recalculates declination nearer the real sunrise time
IF NOT polarday THEN BEGIN
  sundata.time=sp.sunrise
  sp=solarpos(sundata,/LIMB)
  upperL_sunrise=sp.sunrise
  sundata.time=sp.sunset
  sp=solarpos(sundata,/LIMB)
  upperL_sunset=sp.sunset
  ;; real sunrise and sunset
  hour=FIX(upperL_sunrise)
  min=ROUND((upperL_sunrise-hour)*60)/100.
  min=strmid(strtrim(string(min),2),2,2)
  upperL_sunrise=STRTRIM(STRING(hour),2)+':'+min+'h'
  hour=FIX(upperL_sunset)
  min=ROUND((upperL_sunset-hour)*60)/100.
  min=strmid(strtrim(string(min),2),2,2)
  upperL_sunset=STRTRIM(STRING(hour),2)+':'+min+'h'

XYOUTS,.8,.07,'Sunrise: '+upperL_sunrise,charsize=1.0,color=28671,ALIGNMENT=.5,/norm
XYOUTS,.12,.07,'Sunset: '+upperL_sunset,charsize=1.0,color=28671,ALIGNMENT=.5,/norm
ENDIF

end
 ELSE:
  ENDCASE

END

PRO plotsun

COMMON plotsun_block, latitude,longitude,timezone,yearw,julianday,interval,$
    sundata,outfname

;**initial values set to present date
today=BIN_DATE(SYSTIME())
thisyear=FIX(today[0])
thisday=FIX(today[2])
thismonth=FIX(today[1])
thisjulday=JULDAY(thismonth,thisday,thisyear)-JULDAY(1,0,thisyear)

  base = WIDGET_BASE(title='Solar path', /row)
  WIDGET_CONTROL, /MANAGED, base
  b1 = WIDGET_BASE(base, /frame, /column)
  a1 = WIDGET_BASE(b1, /row, /ALIGN_RIGHT)

 help  = widget_button(a1, xsize=50,value = "Help",UNAME="HELP")

  blanks='                           '
 t1 = WIDGET_TEXT(b1, xsize=75, ysize=1,$
    value=blanks+'Solar path: input variables',/ALIGN_CENTER )
 a2 = WIDGET_BASE(b1, /row, /ALIGN_LEFT)

  latitude = cw_field(a2,xsize=15,ysize=1,value='45.93',$
    title='Latitude:   ',/float,UNAME="COMPUTE",/RETURN_EVENTS)

  longitude = cw_field(a2,xsize=15,ysize=1,value='6.87',$
    title='Longitude:               ',/float,UNAME="COMPUTE",/RETURN_EVENTS)

 a3=WIDGET_BASE(b1, /row, /ALIGN_LEFT)

  timezone = cw_field(a3,xsize=15,ysize=1,value='1',$
    title='Time zone:  ',/float,UNAME="COMPUTE",/RETURN_EVENTS)

  interval = cw_field(a3,xsize=15,ysize=1,value='10',$
    title='time interval (minutes): ',/float,UNAME="COMPUTE",/RETURN_EVENTS)

 a4=WIDGET_BASE(b1, /row, /ALIGN_LEFT)

  yearw = cw_field(a4,xsize=15,ysize=1,value=thisyear,$
    title='Year:        ',/float,UNAME="COMPUTE",/RETURN_EVENTS)

  julianday = cw_field(a4,xsize=15,ysize=1,value=thisjulday,$
    title='Day of the Year:       ',/integer,UNAME="COMPUTE",/RETURN_EVENTS)

a5 = WIDGET_BASE(b1, /row, /ALIGN_RIGHT)

  compJD = widget_button(a5, value = "Compute day of the year", uvalue = "computeJD", $
    UNAME="COMPUTEJD",/ALIGN_RIGHT)

  outfname = cw_field(b1,xsize=15,ysize=1,value='solpos.dat',$
    title='Output file name: ',/string,UNAME="COMPUTE",/RETURN_EVENTS)


  compute = widget_button(b1, value = "Compute", uvalue = "COMPUTE", $
    UNAME="COMPUTE",FRAME=5, xsize=150,/ALIGN_CENTER)

  done = widget_button(b1, value = "DONE", uvalue = "DONE", $
    UNAME="DONE",FRAME=5, xsize=150,/ALIGN_CENTER)

widget_control, base, /realize
XMANAGER, 'plotsun', base, /NO_BLOCK


END