test        clear
            set x=0,5,.1
            set y=x*x
            set y[30]=4.5
            set y[2]=2
            limits x y
            box
            con x y
            set sigma=x*0+1
            parabola x y sigma
            ctype red
            con x y_out
            ctype default
            DEL x y sigma x x

            
parabola  3 #Written by C.Austin.  Univ. of Minn. 2005.
            #http://www.astro.umn.edu/~caustin/sm.html
            #this macro fits a parabola using the lsq-2D method
            #define vectors of x values, y values, and error in y values
            #if error in y is unknown, set to unity: "set error=x*0+1"
            #call parabola with 3 arguments: x y y_error
            #this will then output the values of a,b,c,and \chi**2
            #in the variables afit,bfit,cfit,and chisq
            #also will output vector y_out with the fitted yvalues
            set x_in=$1
            set y_in=$2
            set sig_y=$3
            
            #initializing            
            set y_out=x_in*0
            set chisq=x_in*0
            set sum1=0
            set sum2=0
            set sum3=0
            set sum4=0
            set sum5=0
            set sum6=0
            set sum7=0
            set sum8=0
            
            #'math' from numerical recipes...
      	set S_n=sig_y**(-2)
            set S_x=S_n*x_in
            set S_y=S_n*y_in
            set S_xx=S_x*x_in
            set S_xxx=S_xx*x_in
            set S_xxxx=S_xxx*x_in
	      set S_xxy=S_xx*y_in
	      set S_xy=S_x*y_in
            do i=0,dimen(x_in)-1,1{
	      set sum1=$(sum1[0])+$(S_n[$i]) 
	      set sum2=$(sum2[0])+$(S_x[$i])  
	      set sum3=$(sum3[0])+$(S_y[$i]) 
	      set sum4=$(sum4[0])+$(S_xx[$i]) 
	      set sum5=$(sum5[0])+$(S_xxx[$i]) 
	      set sum6=$(sum6[0])+$(S_xxxx[$i]) 
	      set sum7=$(sum7[0])+$(S_xxy[$i]) 
	      set sum8=$(sum8[0])+$(S_xy[$i])}
            set ptemp=$(sum2[0])*$(sum5[0])-$(sum4[0])*$(sum4[0])
            set qtemp=$(sum2[0])*$(sum6[0])-$(sum4[0])*$(sum5[0])
            set rtemp=$(sum4[0])*$(sum6[0])-$(sum5[0])*$(sum5[0])
            set ntemp=$(sum1[0])*$(sum6[0])-$(sum4[0])*$(sum4[0])
            set ttemp=$(sum1[0])*$(sum5[0])-$(sum2[0])*$(sum4[0])
            set stemp=$(sum1[0])*$(sum4[0])-$(sum2[0])*$(sum2[0])
            set d0=$(sum1[0])*$(rtemp[0])-$(sum2[0])*$(qtemp[0])+$(sum4[0])*$(ptemp[0])
            set d1=$(sum3[0])*$(rtemp[0])-$(sum8[0])*$(qtemp[0])+$(sum7[0])*$(ptemp[0])
            set d2=$(sum3[0])*$(qtemp[0])-$(sum8[0])*$(ntemp[0])+$(sum7[0])*$(ttemp[0])
            set d3=$(sum3[0])*$(ptemp[0])-$(sum8[0])*$(ttemp[0])+$(sum7[0])*$(stemp[0])
            set afit=$(d3[0])/$(d0[0])
            set bfit=-$(d2[0])/$(d0[0])
            set cfit=$(d1[0])/$(d0[0])
            
            set y_out=$(cfit[0])+$(bfit[0])*x_in+$(afit[0])*x_in*x_in
            do i=0,dimen(x_in)-1,1{
            set chisq=$(chisq[0])+((y_out[$i]-y_in[$i])/S_n[$i])**2}
            echo fit:   ($(sprintf('%5.3f',$(afit[0]))))x^2+($(sprintf('%5.3f',$(bfit[0]))))*x+$(sprintf('%5.3f',$(cfit[0])))
            echo chisq: $(sprintf('%5.3f',$(chisq[0])))

            DEL sum1 sum2 sum3 sum4 sum5
            DEL sum6 sum7 sum8 S_n S_x
            DEL S_xx S_xxx S_xxxx S_xxy S_xy
            DEL S_y d0 d1 d2 d3
            DEL ntemp ptemp qtemp rtemp stemp
            DEL x_in y_in sig_y ttemp i

	
DEL 5       delete $1
            delete $2
            delete $3
            delete $4
            delete $5