; $Id: wave_dir.pro,v 1.4 2017/01/03 15:27:21 mcnutt Exp $ ;WAVE_DIR: estimate the true direction of a height-time path, assuming uniform radial motion in 3D ; ; $Log: wave_dir.pro,v $ ; Revision 1.4 2017/01/03 15:27:21 mcnutt ; Use select_windows instead of set_plot ; ; Revision 1.3 2008/07/01 17:14:43 lee ; Changed default path for h-t data ; ; Revision 1.2 2008/06/25 22:10:05 sheeley ; made interactive to increase robustness ; ; Revision 1.1 2008/03/20 21:28:57 sheeley ; nbr - copied from adam/caryn on birch ; ;written by Caryn Palatchi with help from Adam Herbst (summer 2007) ; modified by NRS sept. 26, 2007. ;function parameters: A = [delta, tau, t0] = [angle, Dsun/Vr, start time] pro tan_ang, x, a, f, pder f = (x-a[2])*cos(a[0])/(a[1]-(x-a[2])*sin(a[0])) if (N_PARAMS() ge 4) then begin pder = fltarr(N_ELEMENTS(x),3) denom = (a[1]-(x-a[2])*sin(a[0]))^2 pder(*,0) =(x-a[2])*(x-a[2]-a[1]*sin(a[0]))/denom pder(*,1) =-1*(x-a[2])*cos(a[0])/denom pder(*,2) =-1*a[1]*cos(a[0])/denom endif end function wave_dir,vr=vr,delta=delta,range=range list = dialog_pickfile(filter = '*.ht', path = '/sjmaps08/ht_data/', GET_PATH=mypath,/MULTIPLE_FILES) for i=0,n_elements(list)-1 DO BEGIN read_ht_file, list[i], tai, ht, pa if(i EQ 0) then begin ;INITIAL GUESSES AND CONVERSIONS Dsun = 215*6.96e+5; if (data_type(vr) EQ 0) then Vr=260;km/s if (data_type(delta) EQ 0) then delta=45 ; ***********nrs change 6/25/08************ line0: read,'Vr (km/s) =',Vr read, 'delta (deg) =', delta ;******************************************* ;delta=;degrees velocity = (ht[1]-ht[0])/(tai[1]-tai[0]);angular delta=!pi*(delta)/180 tau=Dsun/Vr t0=tai[0]-(ht[0]/velocity) a = double([delta,tau,t0]) vr_initial=vr delta_initial=delta*180/!pi t0_initial=utc2str(tai2utc(t0)) endif if (data_type(range) EQ 0) then range='undefined' ;getting time vs. tan(alpha) from data for different alpha ranges if (range EQ 'undefined') then begin x = double(tai) y = double(tan(!pi*ht/180)) endif if (range EQ 'hi1') then begin cut_hi1=where(ht GT 3.28); -1 if only cor2, first element is 0 if hi1 or hi2 cut_hi2=where(ht GT 23.28);-1 if only hi1 cor2, first element minus 1 is first alpha in hi2 if (cut_hi1[0] EQ -1) then element_start=0 else element_start=cut_hi1[0] if (cut_hi2[0] EQ -1) then element_end=n_elements(ht)-1 else element_end=cut_hi2[0]-1 x = double(tai[element_start:element_end]) y = double(tan(!pi*ht[element_start:element_end]/180)) endif if (range EQ 'hi2') then begin cut_hi2=where(ht GT 18.36) if (cut_hi2[0] EQ -1) then element_start=0 else element_start=cut_hi2[0] element_end=n_elements(ht)-1 x = double(tai[element_start:element_end]) y = double(tan(!pi*ht[element_start:element_end]/180)) endif if (range EQ 'hi1,2') then begin cut_hi1=where(ht GT 3.28) if (cut_hi1[0] EQ -1) then element_start=0 else element_start=cut_hi1[0] element_end=n_elements(ht)-1 x = double(tai[element_start:element_end]) y = double(tan(!pi*ht[element_start:element_end]/180)) endif if (i EQ 0) then time=x if ~(i EQ 0) then time=[time,x] if (i EQ 0) then tan_alpha=y if ~(i EQ 0) then tan_alpha=[tan_alpha,y] endfor ;determining range alpha_min=min(atan(tan_alpha)*180/!pi) alpha_max=max(atan(tan_alpha)*180/!pi) if ((alpha_min GT 3.28) and (alpha_max LT 23.28)) then range='hi1' if (alpha_min GT 18.36) then range='hi2' if ((alpha_min GT 3.28) and (alpha_min LT 18.36) and (alpha_max GT 23.28)) then range='hi1,2' if ((alpha_min LT 3.28) and (alpha_max LT 23.28)) then range='cor2,hi1' if ((alpha_min LT 3.28) and (alpha_max GT 23.28)) then range='all3' ;sort data for curvefit? tta=double(fltarr(n_elements(time),2)) tta[*,0]=time tta[*,1]=tan_alpha time=tta[sort(tta[*,0]),0] tan_alpha=tta[sort(tta[*,0]),1] itmax=1000 tol=1e-7 fit=CURVEFIT(time,tan_alpha,weight,a,sigma,chisq=chisq,/double,function_name='tan_ang',iter=iter,itmax=itmax,status=status,tol=tol,yerror=yerror) ;fit=CURVEFIT(x,y,weight,a,chisq=chisq,/double,function_name='tan_ang') delta=180*a[0]/!pi Vr=Dsun/a[1];km/s t0=utc2str(tai2utc(a[2])) print,list print, 'PA:',pa print, 'range:'+string(range) print, 'number of data points:'+string(n_elements(time)) print, 'status:'+ string(status) print,'iter:'+string(iter) print,'chisq:'+string(chisq) print,'vr_initial:'+string(vr_initial) print,'delta_initial:'+string(delta_initial) print,'t0_initial:'+string(t0_initial) print,'vr:'+string(vr) print,'delta:'+string(delta) print,'t0:'+string(t0) ;print, 'delta error'+ string(sigma[0]) ;print, 'vr error', Dsun/sigma[1] indp=time/3600/24-a[2]/3600/24 dp=180*atan(tan_alpha)/!pi dpfit=180*atan(fit)/!pi ;item0='PA:'+string(pa) item0='PA = '+strmid(strtrim(string(pa),1),0,5)+' deg' ;item1='range:'+' '+ string(range) item1='range = '+strtrim(range,1) item2='number of data points:'+string(n_elements(time)) item3='status:'+ string(status) item4='iter:'+string(iter) ;item5='chisq:'+string(chisq) item5='!4v!3'+'2'+' = '+strmid(strtrim(string(chisq),1),0,13) item6='vr_initial:'+string(vr_initial) item7='delta_initial:'+string(delta_initial) item8='t0_initial:'+string(t0_initial) ;item9='vr:'+string(vr) item9='V!Dr!N = '+strmid(strtrim(string(vr),1),0,5)+' km/s' ;item10='delta:'+string(delta) item10='!4d!3'+ ' = '+strmid(strtrim(string(delta),1),0,4)+' deg' ;item11='t0:'+' '+ string(t0) item11='t!D0!N = '+strmid(strtrim(string(t0),1),0,10)+$ ' '+strmid(strtrim(string(t0),1),11,5)+ ' UT' ;item12='alpha(degrees)' item12='!4a!3'+'(degrees)' item14='days since t!D0!N' item13=string(list) window,0 val = 5*fix(max(dpfit)/5)+10 ;plot,indp,dpfit,xminor=12,xtickinterval=1,ytickinterval=5,xtitle=item14,ytitle=item12 plot,indp,dpfit,xminor=12,xtickinterval=1,ytickinterval=5,xtitle=item14,ytitle=item12,$ charsize=1.5,xrange=[min(indp)-0.5,max(indp)+0.5],xstyle=1,yrange=[0,val],ystyle=1 oplot,indp,dp,psym=4 ;items=[item0,item9,item10,item11,item1,item5,item13] items=[item0,item10,item9,item11,item1] items2=[item0,item10,item9,item11,item1,item5] legend,string(items2),/bottom,/right,spacing=1.2,box=0,charsize=1.5 ;save the results ans='' read, 'want to save(y/n)? ', ans if (ans eq 'y') then begin set_plot,'ps' ans2='' read,'enter filename= ',ans2 device,filename=ans2 ;device,ysize=9,/in device, ysize=7,/in device,xsize=7,/in ;device, yoffset=0.75,/in device, yoffset=2.75,/in ;plot,indp,dpfit,title=ans2,xminor=12,xtickinterval=1,ytickinterval=5,xtitle=item14,ytitle=item12 plot,indp,dpfit,title=ans2,xminor=12,xtickinterval=1,ytickinterval=5,xtitle=item14,$ ytitle=item12, charsize=1.5, xrange=[min(indp)-0.5,max(indp)+0.5],$ xstyle=1,yrange=[0,val],ystyle=1 oplot,indp,dp,psym=4 legend,string(items2),/bottom,/right,spacing=1.2,box=0,charsize=1 device,/close select_windows ;set_plot,'x' ;print the results ans3='' read, 'want to print(y/n)? ', ans3 if (ans3 eq 'y') then begin spawn, 'lp -d hp4-136 '+ string(ans2) endif endif ;**************nrs change 6/25/08****************** ans7='' read,'want to try again (y/n) ',ans7 if (ans7 eq 'y') then goto, line0 ;********************************************* stop wdelete return, a END