*
*
* 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 )'
* 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]