*
*
*  To copy this file to your own machine, either capture the entire
*  text with a mouse, or use the "file/save as" menu to write it out.
*  Make sure you call it example.kumac, and that you write it out as
*  plain text, if possible.  To execute it, just type "exec example" 
*  in PAW.
*  
*
*


Message 
Message 
Message 'Hi - Please read the comments in the kumac.'
Message 
Message 

* This kumac gives some examples of common paw commands.
* 
* It simulates an experiment of measuring the trajectory of 
* an object shot upwards at 100m/s in a vacuum.
* The simulated data points are randomized versions of the ideal
* trajectory.  The simulated points are fit to a parabola, to recover
* the the initial position and velocity, and g.

* The points are randomized with an error that is a constant in time,
* and with an error on the height proportional to the height.


* Set some variables to describe the case (such variables can only be
* set with a kumac)
* When you set a variable, you refer to it as var.  When you use its
* value, you refer to it as [var].
* See help.  Choose the MACRO menu, then SYNTAX, then VARIABLES, etc.

g=-9.8
x0=0
v0=100.
errort=.2

* Set 1 plot per screen
zone

* Book a function, with the ideal height vs. time relationship.
*
* fun1 - the command FUNCTION/FUN1 - creates a 1-d function
*   20 - histogram ID to store the function in.  You can choose any
*        number<9999999(?)
*   [x0]+[v0]*x+[g]*x**2 - the function.  [x0], [v0], and [g] refer to
*                          the variables I set above.
*   10  - Make the histogram 10 bins wide, from x=.5 to x=10.5

fun1 20 [x0]+[v0]*x-abs([g])*x**2 10 .5 10.5

* Label the axis on the plot
atit 'Time (s)' 'Height (m)'

message
wait 'This is the ideal height-time relationship (Hit <Return>)'

* open a new "RZ" file, and write out the histogram containing the function.
* please see the help page in paw!

his/file 1 tmp_example.his ! N
cd //lun1

* This puts histogram 20, which we defined above, into the new file
* tmp_example.his

hrout 20
close 1

* Delete all histograms to prove that it's really saved.
his/del 0

* List all histograms
his/lis

* Write a blank line to the command window.
message
wait 'This is a listing of the currently defined histograms (there are none)'

* Load in the file we just wrote
his/file 1 tmp_example.his 
his/lis

message
wait 'Now we''ve re-loaded the histogram that we had saved to the file'
his/plo 20
close 1

* Delete the RZ file, since we don't really need it. 
* "sh" is short for "shell," which means "execute the rest of the line
* in an operating system shell."
sh rm tmp_example.his



* See if any of the vectors used here exist, and delete them all if so.
* $VEXIST is one of many useful "KUIP" functions.  See the help page
* for kuip/function.

if ($VEXIST(intime).ne.0) then
   vec/del intime
   vec/del iny
   vec/del timeout
   vec/del yout
   vec/del dummy
   vec/del evy
   vec/del evt
endif
 
* Create a vector of errors on the data which will be useful later
vec/crea evy(10) R 
* Create a vector of errors on the time measurement, and fill each of
* the 10 elements with the value of the variable errort.
vec/crea evt(10) R 10*[errort]

* ideal time
vec/crea intime(10) 

* ideal height
vec/crea iny(10)

* simulated "measured" time
vec/crea timeout(10)

* simulated "measured" height
vec/crea yout(10)


* get the ideal vectors from the function we created earlier
get_vec/con 20 iny
get_vec/abs 20 intime

* We need this vector because of the way PAW calls FORTRAN routines.
vec/crea dummy(1) R 1.

* Step through the 10 simulated data points.
do istep=1,10

* Pick a random number to smear the time measurement.  See the macro ranorm
* in this file
 
   exec example#rannorm
* [@] is the value returned by the last macro to be executed.
   dt=[@]
   time=[istep]+[dt]*[errort]
   vec/inp timeout([istep]:[istep]) [time]

* Smear the height measurement (to 10% of itself)
   exec example#rannorm
   dy=[@]
   y=iny([istep])

* SIGMA is a CERN application that is very useful for doing arithmatic in
* PAW.
   errory=$SIGMA(sqrt(([y]*.10)**2+(1.)**2.))
   y=[y]+[dy]*[errory]
   vec/inp yout([istep]:[istep]) [y]
   vec/inp evy([istep]:[istep]) [errory]
   
* this closes the "do istep=1,10" above
enddo


* Book a histogram for the simulated "measured" data.
1dhis 11 'Measured Y vs. T' 10 .5 10.5

* Put the simulated vectors in, as well as errors.
put_vec/con 11 yout
put_vec/err 11 evy

* fit the histogram to a parabola.  See the help page.
his/fit 11 p2

* Note: you could also fit the vectors directly.

message
message 'The parameter corresponding to "g" is A2 in the list'
* The wait command pauses for the user to hit "return"
wait 'This is a fit of a parabola to the simulated data'

* Set the plot symbol to a filled circle.
igset mtyp 20

* overlay the input data plots on the histogram (this isn't stored in the
* histogram, just display on top of it.)
vec/plo iny%intime ! s

* Set the plot symbol back to a dot.
igset mtyp 1
wait 'The dots are the true points'

* Set the window to have one plot above another.
zone 1 2

* Set up a coordinate system for the following commands.
null 0. 12. 1. 300.
* plot the output vectors.
hpl/err timeout yout evt evy 10 26
wait 'Here is another way to plot the data.'

* Set the y axis to be logarithmic
opt logy
null 0. 12. 1. 300.
hpl/err timeout yout evt evy 10 26
wait 'Here are the data again, on a log scale.'
opt liny

* make a 2d function a plot it 2 ways.
fun2 30 [x0]+y*x-abs([g])*x**2 100 .5 10.5 100 90. 110.
his/plo 30 zcol
atit 'Time (s)' 'V0 (m/s)'
message 'These are 2-d plots showing the height-time relationship vs.'
message 'initial velocity'


message
message
message 'Excercise: Run this 20 times, keeping track of each "A2" returned'
message 'by the fit.  Put them into a vector, and vec/plo the vector.'
message 'If the "stat" option is on (see help opt) then the RMS should'
message 'be the error you typically had on "A2".'

* Peter Shanahan - June 5, 1998.


return

*************************************
macro rannorm
*=======================================================================
*     Produces randomly generated numbers according to an normal
*     distribution.  Algorithm by Press, Flannery, Teukolsky & Vetterling.  
*     This macro based on FORTRAN  taken from L. Bellantoni

* This is what a statement label looks like in kuip.
909:

* 
* see "help func" for a list of very useful functions, such as $CALL.
          V1 = 2.0 * $CALL('RNDM(dummy)') - 1.0
          V2 = 2.0 * $CALL('RNDM(dummy)') - 1.0

          R = $SIGMA(([v1])**2+([v2])**2)
          IF ([R].GE.1.0) then
              GOTO 909
          endif
          FAC = $SIGMA(SQRT(-2.0 * LOG([R])/[R]))
          RANORM = [V2]*[FAC]

* You can only return one variable in a PAW macro.
return [ranorm]