;*************************************************
; regress_6.ncl
;
; Concepts illustrated:
;   - Calculating the least squares line fit via 'lspoly_n'           
;   - Calculating the simple linear regression via 'regline_stats'
;   - Drawing a scatter plot and regression and polynomial curves
;*************************************************
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"


begin

;--------------------------------------------------------
; Open Stations
   diri = "./"
   nxy    = numAsciiRow("DAILY.csv")               ; # values
;print(nxy)
;exit   


   xy = asciiread("DAILY.csv",-1, "string")

;printVarSummary(xy)
;print(xy)
;exit
 
   str_xy = str_split_csv(xy,",",3)
   str_xy@_FillValue = -999
;print(str_xy)
;exit

do it=0,1   

if(it .eq. 0) then
   x = tofloat(str_xy(:,4)) ;
   y    = tofloat(str_xy(:,6))   ; Ozone_Megan
   color = "green"
   simulation = "M2.04"
else
   x = tofloat(str_xy(:,5)) ;
   y    = tofloat(str_xy(:,7))   ; Ozone_Megan
   color = "red"
   simulation = "M2.10"
end if


   x@_FillValue = 0
   y@_FillValue = 0
   nxy  = dimsizes(y)
;print(x)
;print(y)
;exit
    
;*************************
; least squares: x,y pairs can be in any order
;**************************
   nc   = 4                            ; 3rd degree polynomial
   coef = lspoly_n(x, y, 1, nc, 0)     ; Least Squares Fit; all weights are set to one 
   print(coef)
   print("======")
;exit
;*************************
; linear regression: x,y pairs can be in any order
;**************************

   rc = regline (x,y)                  ; linear regression
   print(rc)
   print("======")
;exit

;**************************
; Plot lines require that the ascissa be in monotonic order (ascending or descending)
; Note: markers do *not* require any ordering
;**************************

   mono = 1                            ; ascending=1 , descending=-1
   ii   =  dim_pqsort_n(x,mono,0)
;  print(ii)
;exit

   xx   = x(ii)                        ; ascending order 
;print(xx)
   yy   = y(ii)
;print(yy)
;exit
;**************************
; generate the lines
;**************************

   ypoly = coef(0) + coef(1)*xx + coef(2)*xx^2 + coef(3)*xx^3  ; polyline
   yregr = rc*xx + rc@yintercept                               ; regression

;**************************
; PLOT
;**************************

  tplt      = new((/3,nxy/), typeof(x), x@_FillValue)
  tplt(0,:) = yy    
;  tplt(1,:) = ypoly
  tplt(2,:) = yregr


  res                     = True                   ; plot mods desired
  res@gsnFrame            = False     ; Don't advance the frame.
  res@gsnDraw             = False                    ; don't draw plot
  res@gsnMaximize         = True                   ; maximize plot in frame
  res@xyMarkLineModes     = (/"Markers","Lines", "Lines"/)  ; choose which have markers
  res@xyMarkers           = 16                     ; choose type of marker 
  res@xyMarkerColor       = color
  res@xyMonoLineColor     = False
  res@xyLineColors        = (/"green","blue","black"/)
  res@xyMarkerSizeF       = 0.005                  ; Marker size (default 0.01)
  res@xyDashPatterns      = 1                      ; solid line 
  res@xyLineThicknesses   = (/1,5,5/) 
  res@vpWidthF            = 0.75
  res@vpHeightF           = 0.5

  res@trXMinF = 0
  res@trXMaxF = (10) 
  res@trYMaxF = (100)
;  res@trXMaxF = 180
  res@trYMinF = 0
 ; res@trYMaxF = 180
  
  res@tiYAxisString       = "O3"
  res@tiXAxisString       = "HCHO/NO2"
 ; res@tiMainString        =  simulation  ; title
 
type = "x11"
;type = "png"
;type = "pdf"
wks = gsn_open_wks(type,"regress"+simulation)              ; send graphics to PNG file

   plot  = gsn_csm_xy (wks,xx,tplt,res)             ; create plot
;--------------------------------------------------------------
; Add some text strings to the plot, using X,Y axis coordinate
; values to indicate the locations of the strings. The gsn_add_text
; function just attaches the text strings to the plot; it doesn't
; draw them. The text strings won't get drawn until you draw the plot.
; Each text string variable must be unique.
;
   txid = new(2,graphic)                  ; Nine text strings

   txres               = True             
   txres@txFontHeightF = 0.02
   txres@txFontColor = "black"

   txres@txJust        = "CenterLeft"
   txres@txFuncCode    = "~"

   txid(0) = gsn_add_text(wks,plot,"N pt. = "+nxy, 100,30,txres)
   txid(1) = gsn_add_text(wks,plot,"Y = "+sprintf("%5.2f", rc)+"X + "+sprintf("%5.2f", rc@yintercept), 100,15,txres)

   draw(plot)
   frame (wks)

end do
end