function demacs,logt,dem0=dem0,logt0=logt0,norm=norm,slope=slope,$ group=group,igroup=igroup, _extra=e ;+ ;function demacs ; a DEM editor that allows interactively "tweaking" a specified DEM ; ;syntax ; dem=demacs(logt,dem0=dem0,logt0=logt0,norm=norm,slope=slope,$ ; group=group,igroup=igroup, PLOT_KEYWORDS) ; ;parameters ; logt [INPUT; required] log10(Temperature [K]) at which DEM ; is to be determined ; ;keywords ; dem0 [INPUT] initial DEM (just for guidance). ; * if size does not match LOGT0, gets stretched to suit. ; * overrides SLOPE, but not NORM ; * if max(DEM0)<100, assumed to be in log form -- i.e., ; plots are in "linear" form and +vity is not enforced. ; logt0 [INPUT] logT at which DEM0 are defined. if not given, ; assumed to be LOGT. ; norm [INPUT] if set, first element of initial DEM is set to NORM ; slope [INPUT; default=1] if given, generates initial ; DEM=NORM*T0^(SLOPE) ; * if DEM0 is specified, SLOPE is ignored ; group [OUTPUT] long-int array of same size as LOGT containing ; grouping information ; igroup [OUTPUT] array of same size as GROUP, containing index ; of "representative element" of each group ; _extra [INPUT] allows passing defined keywords (e.g., YRANGE) ; to PLOT. ; ;restrictions ; requires X-windows, display, and mouse capability. ; ;history ; vinay kashyap (Feb97) ; corrected log behavior, added _extra to PLOT (VK; Apr97) ; cleaned up log behavior, added keystrokes *,/,v,^,z (VK; Aug98) ; button press status now stored in !MOUSE, not !ERR (VK; Apr09) ;- ; usage nt=n_elements(logt) if nt lt 3 then begin print,'Usage: dem=demacs(logt,dem0=dem0,logt0=logt0,norm=norm,slope=slope,$' print,' group=group,igroup=igroup,PLOT_KEYWORDS)' print,' DEM editor' & return,0. endif ; check keywords ; nt0=n_elements(logt0) ;is LOGT0 defined? if nt0 eq 0 then begin ;nope. lt0=logt & nt0=nt ;use LOGT endif else lt0=[logt0] ;make sure 'tis an array ; if keyword_set(dem0) then d0=[dem0] else begin ;DEM0 as an array if n_elements(slope) eq 0 then slope=1. ;no DEM0, so what slope? d0=10.D^(slope*lt0+10) ;compute DEM0 endelse if n_elements(d0) eq 1 then d0=dblarr(nt0)+d0(0) nd0=n_elements(d0) ;if max(d0) lt 100 then d0=10.^(d0) ;in log form? if max(d0) lt 100 then ylog=0 else ylog=1 ;in log form? ; NOTE: ylog=0 says that the *input* is in LOG form! if nd0 ne nt0 then begin t0=lt0(0) & t1=lt0(nt0-1) & ti=findgen(nd0)*(t1-t0)/(nd0-1.)+t0 d0=interpol(d0,ti,lt0) ;stretch to fit endif ; if keyword_set(norm) then d0=d0*norm/d0(0) ;input normalization ; ; groups of points if n_elements(group) eq nt then grp=group else grp=lindgen(nt) ; representative point in each group if n_elements(igroup) eq nt then igrp=igroup else igrp=grp ; initialize mbt=0 ;buttons pressed kbrd=0 ;input from keyboard go_on=1 ;active mouse act='?' ;action to take factor=1.1 ;multiplicative "nudge" factor togint=1 ;spline/poly/linear interpolation ;dem=interpol(d0,lt0,logt)>1 ;output dem=(spline(lt0,d0,logt,1.)>min(d0))0 ; hop left 'l': ip=(ip+1)<(nt-1) ; hop right ' ': ip=(ip+1)<(nt-1) ; hop right 'b': ip=(ip-4)>0 ; skip left 'w': ip=(ip+4)<(nt-1) ; skip right 'B': ip=(ip-16)>0 ; jump left 'W': ip=(ip+16)<(nt-1) ; jump right '0': ip=0 ; move to 1st element '$': ip=nt-1 ; move to last element 'j': begin ; divide if ylog eq 1 then dd(igrp(ip))=dd(igrp(ip))/factor else $ dd(igrp(ip))=dd(igrp(ip))-alog10(factor) recomp=1 end 'J': begin ; divide if ylog eq 1 then dd(igrp(ip))=dd(igrp(ip))/2 else $ dd(igrp(ip))=dd(igrp(ip))-alog10(2.) recomp=1 end '-': begin ; divide if ylog eq 1 then dd(igrp(ip))=dd(igrp(ip))/10 else $ dd(igrp(ip))=dd(igrp(ip))-1. recomp=1 end 'k': begin ; multiply if ylog eq 1 then dd(igrp(ip))=dd(igrp(ip))*factor else $ dd(igrp(ip))=dd(igrp(ip))+alog10(factor) recomp=1 end 'K': begin ; multiply if ylog eq 1 then dd(igrp(ip))=dd(igrp(ip))*2 else $ dd(igrp(ip))=dd(igrp(ip))+alog10(2.) recomp=1 end '+': begin ; multiply if ylog eq 1 then dd(igrp(ip))=dd(igrp(ip))*10 else $ dd(igrp(ip))=dd(igrp(ip))+1. recomp=1 end '*': begin ; multiply whole group og=where(grp eq grp(ip),mog) if mog eq 0 then message,'bug!' if ylog eq 1 then dd(og)=dd(og)*factor else $ dd(og)=dd(og)+alog10(factor) olddd=dem & dem=dd end '/': begin ; divide whole group og=where(grp eq grp(ip),mog) if mog eq 0 then message,'bug!' if ylog eq 1 then dd(og)=dd(og)/factor else $ dd(og)=dd(og)-alog10(factor) olddd=dem & dem=dd end 'v': begin ; divide whole group og=where(grp eq grp(ip),mog) & tmp=0. if mog eq 0 then message,'bug!' while tmp le 0 do $ read,prompt='divide this group by what factor? ',tmp if ylog eq 1 then dd(og)=dd(og)/tmp else $ dd(og)=dd(og)-alog10(tmp) olddd=dem & dem=dd end '^': begin ; multiply whole group og=where(grp eq grp(ip),mog) & tmp=0. if mog eq 0 then message,'bug!' while tmp le 0 do $ read,prompt='multiply this group by what factor? ',tmp if ylog eq 1 then dd(og)=dd(og)*tmp else $ dd(og)=dd(og)+alog10(tmp) olddd=dem & dem=dd end 'f': begin ; x by... read,prompt='nudge factor: ',factor if factor eq 0 then factor=1.1 if factor lt 1 then factor=abs(1./factor) end 'F': begin ; ask what value tmp=0. while tmp le 0 do $ read,prompt='DEM(logT='+strtrim(logt(igrp(ip)),2)+') =? ',tmp dd(igrp(ip))=tmp recomp=1 end 'g': begin ; fix value dd(ip)=y0 recomp=1 end 'n': begin ; normalization tmp=0. while tmp le 0 do $ read,prompt='set normalization @ logT='+strtrim(logt(ip),2)+': ',tmp if ylog eq 1 then dd=dd*tmp/dd(ip) else dd=dd+tmp-dd(ip) recomp=1 end '(': begin ; start grouping ig0=ip & act2='' while act2 ne ')' do begin ;(group act2=get_kbrd(1) case act2 of 'h': ip=(ip-1)>0 ; move left 'l': ip=(ip+1)<(nt-1) ; move right ' ': ip=(ip+1)<(nt-1) ; move right( else: print,'h: left; l: right; ): exit' endcase oplot,[logt(ip)],[dd(ip)],psym=6 & wait,0.02 oplot,[logt(ip)],[dd(ip)],psym=6,col=0 oplot,[logt(ip)],[dd(ip)],psym=7 endwhile ;end group) act='>' & ig1=ip end '{': begin ; start ungrouping jg0=ip & act2='' while act2 ne '}' do begin ;(ungroup act2=get_kbrd(1) case act2 of 'h': ip=(ip-1)>0 ; move left 'l': ip=(ip+1)<(nt-1) ; move right ' ': ip=(ip+1)<(nt-1) ; move right{ else: print,'h: left; l: right; }: exit' endcase oplot,[logt(ip)],[dd(ip)],psym=6 & wait,0.02 oplot,[logt(ip)],[dd(ip)],psym=6,col=0 oplot,[logt(ip)],[dd(ip)],psym=7 endwhile ;end ungroup) act=']' & jg1=ip end 'u': begin ; undo last change dem=olddd & olddd=dd ;recomp=1 end 'U': begin ; ungroup just one point jg0=ip & jg1=ip & act=']' end 't': begin ; toggle spline/poly/lin print,'1:cubic spline, 2:polynomial, 3:linear' act2='' & act2=get_kbrd(1) if act2 eq ' ' or act2 eq '' then togint=togint+1 else togint=fix(act2) recomp=1 end '.': gg=ip ; representative point '?': print,hilfe ; help 'x': kbrd=0 ; cursor 'q': go_on=0 ; quit 'z': stop,'halting' ; stop '>': ; same as group ']': ; same as ungroup else: print,'?: help; x: cursor; q: quit' endcase ;action} oplot,[logt(ip)],[dem(ip)],psym=4 if act eq '>' then begin ;group if ig1 gt ig0 then begin grp(ig0:ig1)=grp(ig0) igrp(ig0:ig1)=igrp(ig0) oplot,logt(ig0:ig1),dem(ig0:ig1),psym=1 endif endif if act eq '.' or act eq 'g' then begin ;representative point oo=where(grp eq grp(ip)) igrp(oo)=gg oplot,logt(ig0:ig1),dem(ig0:ig1),psym=4,col=0 oplot,logt(ig0:ig1),dem(ig0:ig1),psym=1 oplot,[logt(ip)],[dem(ip)],psym=4 endif if act eq ']' then begin ;ungroup if jg1 ge jg0 then begin grp(jg0:jg1)=lindgen(jg1-jg0+1)+jg0 igrp(jg0:jg1)=grp(jg0:jg1) oplot,logt(jg0:jg1),dem(jg0:jg1),psym=1,col=0 oplot,logt(jg0:jg1),dem(jg0:jg1),psym=4,col=0 oplot,logt(jg0:jg1),dem(jg0:jg1) if jg1 lt nt-1 and jg0 gt 1 then begin if grp(jg0-1) eq grp(jg1+1) then begin o1=where(grp eq grp(jg0-1) and lindgen(nt) lt jg0) o2=where(grp eq grp(jg1+1) and lindgen(nt) gt jg1) grp(o2)=jg1+1 if igrp(jg0-1) ge jg0 then igrp(o1)=jg0-1 if igrp(jg1+1) le jg1 then igrp(o2)=jg1+1 endif endif endif endif ilogt=igrp(uniq(igrp,sort(igrp))) if n_elements(ilogt) lt 3 then recomp=0 if recomp eq 1 then begin ;(recompute lt0=logt(ilogt) & d0=dd(ilogt) if total(abs(dem(ilogt)-dd(ilogt)))/total(dd(ilogt)) gt 1e-7 then olddd=dem if togint gt 4 then togint=1 case togint of 1: begin ;cubic spline ;dem=spline(lt0,d0,logt,0.1) dem=spline(lt0,d0,logt,1.) end 2: begin ;polynomial if ylog eq 1 then dem=10.D^(spline(lt0,alog10(d0),logt,5)) else $ dem=spline(lt0,d0,logt,5) end 3: begin ;linear dem=interpol(d0,lt0,logt) end else: begin ;spline interpolate y2=spl_init(lt0,d0,/double) dem=spl_interp(lt0,d0,y2,logt,/double) end endcase oo=where(dem gt 0,moo) if moo gt 0 then begin if ylog eq 1 then dem=dem>(min(dem(oo))) dd=dem endif else begin dd=olddd & dem=dd endelse endif ;recomputed) endwhile ;go_on=1} group=grp & igroup=igrp & return,dem ;outputs end