;The procedure phitheta_fit_wc.pro takes Willow Creek rotation data and fits ;it with a function for use in combineddayflux. Phi is the wind direction and ;theta is the tilt angle. Both are input in radians. This code should work for ;rot files from dayflux_wc.pro concatenated with create_annual_wc.pro. ; ;Use is as follows: ; phitheta_fit_wc,plt,yyyy,medfilt,arr,newamin,newomega,co ;where ; | plt = 1 to plot the data and resulting fit to a ; | postscript file: ; | /data/davis/cheas/wcreek/data/rot/rot_son.ps ; | where 'son' is a string described below ; | yyyy = 4-digit year ; input | medfilt = 1 to use median filtered angles for fit ; | 0 to not use median filtered angles for fit ; | son = a string denoting the level and time period ; | for the sonic (e.g., 'son6_5') ; | arr = input array (wcyyyy.rot without header) ; | Format expected is that given by ; ; | newamin = magnitude of fitted sine function ; output | newomega = phase of fitted sine function ; | co = coefficient matrix from polynomial fit ; ; The fitted function is: ;theta_fitted=aminnew*sin(phi-omeganew)+co(3)*phi^3+co(2)*phi^2+co(1)*phi+co(0) ; ;The program works basically as follows: The operator is asked to look at a ;plot of the data and to give an appoximate phase. Then the maginitude of an ;optimized (error is minimized) sine function is fitted. Then the phase that ;produces the least error is found. Then the sine magnitude is recalculated ;and the phase refound. The program continues to find magnitudes and phases ;until the difference between the error from subsequent steps is sufficiently ;small. ; ;Written by BWBerger 7/19/99 based on phitheta_tst.pro. ;This program allows the phase of the fitted sine function to be optimized. ;Modified on 8/30/99 by BWBerger to add norm and medfilt input flags. ;modified on 10/4/99 by weiguo wang to have it suitable for Wcreek. pro phitheta_fit_lc_highorder,plt,yyyy,medfilt,arr,newamin,newomega,co yyyy=strtrim(string(yyyy),2) if medfilt eq 0 then begin son=yyyy+'_lc_mean' subt='_mean' endif else begin son=yyyy+'_lc_median' subt='_median' endelse ;extract phi and theta from array based on value of medfilt (radians) if medfilt eq 0 then begin ;don't use median filtered angles phi=arr(7,*) theta=arr(6,*) endif else begin ;use median filtered angles phi=arr(9,*) theta=arr(8,*) endelse ;First find where phi and theta are both good goodi=where((phi ne -999.) or (theta ne -999.), count) ;Then create variables containing the good data if (count gt 0) then begin phig=phi(goodi) thetag=theta(goodi) ;Also rough screen the data - get indexes when phi<=30 degrees goodii=where(abs(thetag) le 30.*!pi/180.,count) if count gt 0 then begin phig1=phig(goodii) thetag1=thetag(goodii) ENDIF goodii=where( thetag ge -4.*!pi/180.,count) if count gt 0 then begin phig=phig(goodii) thetag=thetag(goodii) endif ;Define good phi and theta as double precision thetag=double(thetag) phig=double(phig) n=float(len(phig)) ;Initialize amin, omega and newomega amin=0. amin=double(amin) omega=0. omega=double(omega) newomega=0. newomega=double(newomega) ;Plot data and have operator select starting phase value plot,phig1*180./!pi,thetag1*180./!pi,psym=3,$ xrange=[-180,180],xstyle=1,$ yrange=[-30,30],ystyle=1,$ xtitle='!7u!3 (deg)',ytitle='!7h!3 (deg)',$ title='LOST CREEK '+yyyy+subt refline,1,1,0.,0 print,'Enter +/- degrees phase shift:' print,'Note: ~180 degree phase shift is not uncommon.' print,'Positive shifts sine to right, Negative to left.' read,omega ;convert to radians omega=omega*!pi/180. ;initialize phase vector used to find optimal phase and corresponding error ;make phase span +/- 20 degrees of omega. 0.02 deg resolution phase=(findgen(2001)/50.-20)*!pi/180.+omega pherror=fltarr(2001) ;initialize an angle vector for plotting purposes angle=(findgen(361)-180)*!pi/180. ;loop while error is unacceptable or number of iterations is low stopflag=0 i=0 while (stopflag eq 0) do begin ;calculate magnitude of sine for min error (use minimization equation) amin=total(thetag*sin(phig-omega))/total(sin(phig-omega)^2) ;calculate phase for min error by checking error for nearby phase values for j=0,len(phase)-1 do begin pherror(j)=total((thetag-amin*sin(phig-phase(j)))^2) endfor ;find index of min pherror and then phase at min pherror ;set newomega to optimized phase gpherror=where(pherror eq min(pherror),cpherror) if cpherror ne 0 then begin bestindex=ceil((len(gpherror)-1)/2.) ;centers index if multiple mins. newomega=phase(gpherror(bestindex)) endif else begin print,'Problems finding phase.' goto, bottom endelse ;Re-find amin based on optimized phase (newomega) newamin=total(thetag*sin(phig-newomega(0)))/total(sin(phig-omega(0))^2) ;Now find polynomial fit to residual (theta data minus sine fit) ;I don't understand why newomega ends up being an array, or why ;its subtraction with phig returns a single value unless (0) is used. thetag_sine=newamin*sin(phig-newomega(0)) co=poly_fit(phig,thetag-thetag_sine,8,/double) print,'co=',co thetag_resid=co(8)*phig^8$ +co(7)*phig^7+co(6)*phig^6+co(5)*phig^5+co(4)*phig^4$ +co(3)*phig^3+co(2)*phig^2+co(1)*phig+co(0) thetag_tot=thetag_sine+thetag_resid ;Get the error for the total fitted function error_tot=sqrt(total((thetag-thetag_tot)^2)/n) ;for plotting, calculate the sine function, the residual poly function ;and the total function theta_sine=newamin*sin(angle-newomega(0)) theta_resid=co(8)*angle^8$ +co(7)*angle^7+co(6)*angle^6+co(5)*angle^5+co(4)*angle^4$ +co(3)*angle^3+co(2)*angle^2+co(1)*angle+co(0) theta_tot=theta_sine+theta_resid ;Print out a bunch of calculated values print,'Initial amin=' print,amin print,'Optimized amin=' print,newamin print,'Initial omega=' print,omega*180./!pi,' degrees' print,'Optimized omega=' print,newomega*180./!pi,' degrees' print,'co*180/pi=' print,co*180./!pi if i ne 0 then begin print,'Old total error=' print,old_error_tot endif print,'Total error=' print,error_tot ;Get strings of parameters for plotting label purposes a=strtrim(string(newamin),2) a_deg=strtrim(string(newamin*180./!pi),2) ph=strtrim(string(newomega),2) ph_deg=strtrim(string(newomega*180./!pi),2) ;higher order fitting c8=strtrim(string(co(8)),2) c7=strtrim(string(co(7)),2) c6=strtrim(string(co(6)),2) c5=strtrim(string(co(5)),2) c4=strtrim(string(co(4)),2) c3=strtrim(string(co(3)),2) c2=strtrim(string(co(2)),2) c1=strtrim(string(co(1)),2) c0=strtrim(string(co(0)),2) c0_deg=strtrim(string(co(0)*180./!pi),2) avgdev=strtrim(string(error_tot),2) avgdev_deg=strtrim(string(error_tot*180./!pi),2) ; ns=strtrim(string(long(n))) ns=strtrim(string(long(n/2.0))); in the wc, a sample every half an hour ;Plot the total fitted function on top of the data plot,phig1*180./!pi,thetag1*180./!pi,psym=3,xticks=8,$ xtickname=['-180','-135','-90','-45','0','45','90','135','180'],$ xrange=[-180,180],xstyle=1,$ yrange=[-30,30],ystyle=1,$ xtitle='!7u!3 (deg)',ytitle='!7h!3 (deg)',$ title='LOST CREEK '+yyyy+subt oplot,angle*180./!pi,(theta_tot)*180./!pi,color=7600000,thick=2 refline,1,1,0.,0 xyouts,.12,.9,'!7h = !3'+a+'*SIN[!7u - !3('+ph+')]+',/normal,charsize=1.2 ; xyouts,.12,.85,'('+c3+')*!7u!3!E3!N+('+c2+')*!7u!3!E2!N+('+c1+')*!7u!3+('+c0+') radians',/normal,charsize=1.1 xyouts,.12,.85,'('+c8+')*!7u!3!E8!N+('+c7+')*!7u!3!E7!N+('+c6+')*!7u!3!E6!N+',/normal,charsize=1.1 xyouts,.12,.80,'('+c5+')*!7u!3!E5!N+('+c4+')*!7u!3!E4!N+',/normal,charsize=1.1 xyouts,.12,.75,'('+c3+')*!7u!3!E3!N+('+c2+')*!7u!3!E2!N+('+c1+')*!7u!3+('+c0+') radians',/normal,charsize=1.1 xyouts,.12,.20,'sine magnitude='+a+' rad ='+a_deg+' deg',/normal xyouts,.12,.16,'phase='+ph+' rad ='+ph_deg+' deg',/normal xyouts,.12,.12,'offset='+c0+' rad ='+c0_deg+' deg',/normal xyouts,.62,.20,ns+' hours of data',/normal xyouts,.62,.16,'avg. dev. from function=',/normal xyouts,.62,.12,avgdev+' rad = '+avgdev_deg+' deg',/normal ;See if change in error is small enough to exit this loop. ;Also exit if number of iterations is too large. if i ne 0 then begin print,'old_error_tot - error_tot =' print,abs(old_error_tot-error_tot) print,'0.001*error_tot =' print,0.001*error_tot if (abs(old_error_tot-error_tot) lt 0.001*error_tot) or (i ge 20) then begin stopflag=1 endif endif print,'**************************' i=i+1 old_error_tot=error_tot ;reset old_error to new value omega=newomega ;resent starting phase to new value endwhile endif else begin print,'No good data' endelse if plt eq 1 then begin set_plot,'ps' ; fname='/data/davis/cheas/wcreek/data/rot/rot'+son+'.ps' ; fname='/davis/s1/cheas/lcreek/data/rot/rot'+son+'.ps' ;fname='/davis/s2/lcreek/data/rot/rot'+son+'highorder.ps' fname='/eddy/s4/lcreek/data/rot/rot'+son+'highorder.ps';change directory on 11/3/2003 device,filename=fname,/landscape plot,phig1*180./!pi,thetag1*180./!pi,psym=3,xticks=8,$ xtickname=['-180','-135','-90','-45','0','45','90','135','180'],$ xrange=[-180,180],xstyle=1,$ yrange=[-30,30],ystyle=1,$ xtitle='!7u!3 (deg)',ytitle='!7h!3 (deg)',$ title='LOST CREEK '+yyyy+subt oplot,angle*180./!pi,(theta_tot)*180./!pi,thick=2 refline,1,1,0.,0 xyouts,.12,.9,'!7h = !3'+a+'*SIN[!7u - !3('+ph+')]+',/normal,charsize=1.2 ; xyouts,.12,.85,'('+c3+')*!7u!3!E3!N+('+c2+')*!7u!3!E2!N+('+c1+')*!7u!3+('+c0+') radians',/normal,charsize=1.1 xyouts,.12,.85,'('+c8+')*!7u!3!E8!N+('+c7+')*!7u!3!E7!N+('+c6+')*!7u!3!E6!N+',/normal,charsize=1.1 xyouts,.12,.80,'('+c5+')*!7u!3!E5!N+('+c4+')*!7u!3!E4!N+',/normal,charsize=1.1 xyouts,.12,.75,'('+c3+')*!7u!3!E3!N+('+c2+')*!7u!3!E2!N+('+c1+')*!7u!3+('+c0+') radians',/normal,charsize=1.1 xyouts,.12,.20,'sine magnitude='+a+' rad ='+a_deg+' deg',/normal xyouts,.12,.16,'phase='+ph+' rad ='+ph_deg+' deg',/normal xyouts,.12,.12,'offset='+c0+' rad ='+c0_deg+' deg',/normal xyouts,.62,.20,ns+' hours of data',/normal xyouts,.62,.16,'avg. dev. from function=',/normal xyouts,.62,.12,avgdev+' rad = '+avgdev_deg+' deg',/normal device,/close set_plot,'x' endif cplot bottom:return end