; solutions for exercise 2 // phy 535 // fall 2011
;
; instructions:
;	at the IDL prompt, type
;
;	IDL> .run soln_exercise2
;
;	the proposal will stop a couple of times so that you (and i)
;	can look at the results. when you are done looking at the 
;	plots, type
;
;	IDL> .continue
;
;	you'll have to continue twice to run all the way through
;	the code.
;
; det // 22 oct 2011
;

; open the windows
window,0
window,1
window,2
window,3

; set up some stuff
counter = 1 ; usually IDL starts with zero, but we need to do this for the plot below
meanresid = make_array(15000) ; make it bigger than we need; 11^4 is 14641.
solution = make_array(11,11,11,11) ; 4-space because we have a four-dimensional trial variables

; read the data
readcol,'exercise2.data.dat',time,flux
wset,0
plot,time,flux,psym=4 ; what does the data look like?

; find a guess solution for y=f(x)
; here are the initial guess constants (by eye)
amp = 3.0
phase = 0.0
omega = 13.5
slope = 6./300.
; now create the guessed function and overplot
xx=findgen(280)
yy = (amp * (sin(xx/omega)+phase)) + (slope*xx)
oplot,xx,yy

; now make a grid over variables of interest
; this loop takes a long time (~5 min?) but its pretty cool to watch the plots
for i=0,10 do begin
   ampt = amp+(((-0.5)+(float(i)/10.))*amp)
   for j=0,10 do begin
      phaset = ((-0.5)+(float(j)/10.))+phase
      for k=0,10 do begin
         omegat = omega+(((-0.5)+(float(k)/10.))*omega)
         for m=0,10 do begin
            slopet = ((-0.05)+(float(m)/100.))+slope

; produce the new function (stepped over the grid) and plot the test function over the data
            yytest = (ampt * (sin(time/omegat)+phaset)) + (slopet*time)
            wset,0
            plot,time,flux,psym=4,symsize=1.5,title='data and possible fits'
            oplot,time,yytest,psym=-4

; calculate the residuals and plot those too
            resid = flux-yytest
            wset,1
            plot,time,resid,psym=4,yrange=[-8,8],title='residuals'
            meanresid[counter] = mean(abs(resid))
            solution[i,j,k,m] = meanresid[counter]
            wset,2
            plot,findgen(counter),meanresid[0:counter],psym=3,xrange=[0,15000],yrange=[0,8],title='solutions',ytitle='mean residual'
            counter = counter + 1
         endfor ; close loop over m
      endfor ; close loop over k
   endfor ; close loop over j
endfor ; close loop over i

stop ; take a look at the plots

; now we can plot the results
; this is basically minimizing over the variables

; amplitude
window,0
plot,[0],/nodata,xrange=[1,5],yrange=[0,2],title='amplitude',ytitle='mean residual'
for i=0,10 do $
   plots,amp+(((-0.5)+(float(i)/10.))*amp),min(solution[i,*,*,*]),psym=4

; phase 
window,1
plot,[0],/nodata,xrange=[-0.6,0.6],yrange=[0,2],title='phase',ytitle='mean residual'
for j=0,10 do $
   plots,((-0.5)+(float(j)/10.))+phase,min(solution[*,j,*,*]),psym=4

; omega
window,2
plot,[0],/nodata,xrange=[5,20],yrange=[0,2],title='omega',ytitle='mean residual'
for k=0,10 do $
   plots,omega+(((-0.5)+(float(k)/10.))*omega),min(solution[*,*,k,*]),psym=4

; slope
window,3
plot,[0],/nodata,xrange=[-0.05,0.07],yrange=[0,6],title='slope',ytitle='mean residual'
for m=0,10 do $
   plots,((-0.05)+(float(m)/100.))+slope,min(solution[*,*,*,m]),psym=4

stop ; take a look at the plots

; now plot the best-fit solution 
; having run this before, i know that the best-fit solution is
; $A=2.4$, 
; phase of zero, 
; omega of~13.5, and
; slope of~0.02.
wset,0
plot,time,flux,psym=4 ; 

; best solution
amp = 2.4
phase = 0.0
omega = 13.5
slope = 0.02
; now create the guessed function and overplot
xx=findgen(280)
yy = (amp * (sin(xx/omega)+phase)) + (slope*xx)
oplot,xx,yy

; now calculate the values at the time points in the data
yy = (amp * (sin(time/omega)+phase)) + (slope*time)
; and calculate the residuals here
resid = flux-yy

; and plot the residuals
wset,1
plot,time,resid,psym=4,yrange=[-8,8]

; and a histogram of the residuals
wset,2
result=histogram(resid,binsize=0.5,locations=loc,min=-2)
plot,loc,result,psym=10

end