SULI Log 2003

2003-07-30T17:28:00.000-07:00

Professor Thomas M. Cover

I e-mailed this guy after receiving an e-mail from Chun-Xi Wang regarding Principal Component Analysis.

Some Ideas To Think About
The FFT of the BPM matrix is the key. We can quickly isolate the betatron frequency, so why do we use SVD? Because we don’t care about the frequency—it is only a tool. What yields the orbits is the amplitude. I am exploring SVD because the current MIA model does not SVD enhance data. I would like to explore how much error this causes relative to the “actual” case and the SVD-enhanced case.

We are also now thinking about phase. We have a “toy” (as Dr. Yan called it) to model noisy BPMs. We can explore the betatron phase by tan^-1(sinlike/coslike), or tan^-1(B/A). Note, this is different from the BPM phase that I inserted into my ideal BPM data matrix.

Things to play with in my “toy”: Right now the period is perfect for periodic motion. What if I change the period? What is the significance of the BPM phase and the angular frequency? I might as well change my variable names to better see this in my equations.

What about the eigenvectors in V? What do those mean? What do they mean when they are chopped off?

2003-07-30T11:12:00.000-07:00

INTERESTING! See "explorations using SVD" file that i made!

2003-07-29T14:21:00.000-07:00

Questions:
1. Am I doing the right thing by "chopping" the last mode? Or should I just set it equal to zero?

2003-07-28T15:56:00.000-07:00

Preliminary results (using findBetaAuto):
Bxx=2.3153
Bxy=2.3162
Byx=2.7687
Byy=2.3180

2003-07-28T15:03:00.001-07:00

function [trim]=findTrim(fxx, xbuferXs)
% function [trim]=findTrim(cfxx)
% July 28, 2003 by Philip Tanedo
% originally from findBetaAuto
% made into a separate function
% ...for use with findBetaAutoSC

tempPeak=0;
tempIndex=0;
i=0;
[tempIndex, tempPeak]=findPeak(fxx);
peakVec(i+1)=abs(tempPeak);
for(i=1:15)
temp=xbuferXs;
for(j=1:i)
temp=BotRowChop(temp); % fine-tuning the frequency
end
fxx=fft(temp);
fxx=TopRowChop(fxx); % chop closed orbit
[tempIndex, tempPeak]=findPeak(fxx);
peakVec(i+1)=abs(tempPeak);
end
[tempIndex, tempPeak]=findPeak(peakVec);
trim=tempIndex-1; % because of note above: to accoun for zeroth trim, index is i+1
disp(['TRIM = ', num2str(trim) ]);

2003-07-28T15:03:00.000-07:00

function [beta]=findBetaAuto(xbuferXs, label)
% function [beta]=findBetaAuto(xbuferXs, label)
% July 28, 2003 by Philip Tanedo
% MODIFIED version of findBeta
% Automates as many processes as possible
% Returns betatron oscillation for fxx
% Preconditions:
% ...flipBuf.mat has been loaded, 'tools' added to path
% Procedure:
% 0. fxx=fft(XbuferXs)
% 1. kill closed orbit
% 2. plot all for user
% 3. prompt user for reasonable data
% 4. close all
% 5. plot representation of frequency increments
% 6. prompt user for reasonable size
% 7. adjust size
% 8. isolate peak frequency
% 9. return peak frequency as betatron frequency
% NOTE: can also be used w/ noisy data to extract
% ...noise frequency when noise is bigger than betatron

% 0
fxx = fft(xbuferXs);

% 1
fxx=TopRowChop(fxx);

% 2
% graphAllFC(fxx, label);

% 3
disp('Please select a reasonable plot to use as key');
choice=input('Index: ');

% 4
close all;

% STEPS 5 and 6 from original findBeta function are commented out
% see below for modification (automate)

% % 5
% i=0;
% figure;
% plot(abs(fxx(:,choice)));
% legend(num2str(choice), 0);
% xlabel('Frequency')
% ylabel('Magnitude (complex)')
% title(['BPM frequency content for BPM#s: ' num2str(choice) ' TRIM VALUE: ' num2str(i)])
%
% for(i=1:15)
%
% temp=xbuferXs;
% for(j=1:i)
% temp=BotRowChop(temp); % fine-tuning the frequency
% end
%
% fxx=fft(temp);
% fxx=TopRowChop(fxx); % chop closed orbit
%
% figure
%
% plot(abs(fxx(:,choice)));
% legend(num2str(choice), 0);
% xlabel('Frequency')
% ylabel('Magnitude (complex)')
% title(['BPM frequency content for BPM#s: ' num2str(choice) ' TRIM VALUE: ' num2str(i)])
% end
%
% % 6 -- note: I could automate this step, as well
% disp('Please select a reasonable plot to use as key');
% trim=input('Trim Value: ');

% 5 and 6 automated

trim=findTrim(fxx, xbuferXs);

% 7
close all;

temp=xbuferXs;
for(j=1:trim)
temp=BotRowChop(temp); % fine-tuning the frequency
end
fxx=real(fft(temp));
fxx=TopRowChop(fxx); % chop closed orbit

% 8
% [pindex, pvalue]=findPeak(fxx);

% but now what does this mean at all?
% finding the right frequency should be trivial

% DESIGN NOTE
% I will design this function to find the complex frequency
% I think this should be good enough to find the sin/cos freq
% Otherwise, I can always modify new functions

% Below code is cited from
% "Spectral Analysis using the FFT" by Bret Ninness
% http://www.ee.newcastle.edu.au/brett/elec2400/matlab4.pdf
% Page 2

% Start citation
temp=size(fxx);
ws=2*pi/temp(1);
wnorm=-pi:ws:pi;
wnorm=wnorm(1:temp(1));
% flip's comment: this is a nifty idea!
% I should have incorporated something like that last line in my
% ... "chop" functions
% End

%fsxx: fft shift of fxx
fsxx=fftshift(fxx);
[spindex, spvalue]=findPeak(fsxx);
beta=abs(wnorm(spindex));

2003-07-28T11:51:00.000-07:00

Things to remember: for the discrete fourier transform, there is a max resolution!

2003-07-28T10:47:00.000-07:00

Writing FFT script. A little bit tedious. How do I convert from FFT frequency to "real" frequency?

2003-07-25T17:45:00.000-07:00

Introduction

In May of 2003, physicists and engineers at the Stanford Linear Accelerator Center were able to achieve record-breaking, peak luminosity in the PEP-II positron/electron storage rings. (Woods) This unprecedented improvement was due in part to a statistical technique called Model Independent Analysis (MIA), which evaluates data without a priori information about the underlying physical apparatus. In order to understand the mechanism by which MIA is able to analyze PEP-II beam data, however, it is useful to first be able to model a simplified particle accelerator in terms of linear optics.

In such a representation, we consider the trace of a particle’s transverse position (x and y coordinates as the particle travels in the z-direction) in 4-dimensional phase space (x and y position along with respective momenta). By representing each lattice element, such as a focusing-defocusing quadrupole magnet, by a matrix transformation acting on a phase space vector, we can extract turn-by-turn maps for the particle’s phase-space position. One can decompose these turn-by-turn maps, or local transfer matrices, into simple rotations in normalized, decoupled phase space, which makes the beam dynamics much simpler mathematically.

Unfortunately, in the real machine we are unable to directly quantify particles’ momenta, and are limited to the data acquired from beam position monitors (BPMs) to determine the center of the electron beam. With what appears to be a mathematical sleight of hand, however, one may still extract the local transfer matrix from this information. BPM data along a given axis (x or y) may be represented as an n x m matrix where n is the number of consecutive-turn measurements and m is the number of BPMs. By exciting the beam’s two eigenplanes and extracting sine-like and cosine-like orbits for each excitation, one collects four linearly-independent orbits which completely determine the linear optics of the system. For example, the data from the excitation of the beam’s first eigenplane results in two data matrices, one representing the x position of the beam at each x-sensitive BPM and the other representing the corresponding y position at each y-sensitive BPM. (In PEP-II not every BPM records both x and y data.) As eigenplane one lies closer to the x-axis of the BPMs, it is a much clearer source of information about the sine-like and cosine-like frequency of betatron oscillations.

BPM noise is removed from the data through the use of FFT and SVD, which also provide a way to identify noisy BPMs and even extract dispersion from synchrotron modes. To fit data, one assumes the longitudinal position of the beam’s magnets are fixed and known. BPM errors and the precise strength of each magnet are then treated as variables and a numerical algorithm is iterated to fit the clean BPM data with those variables. The machine’s ideal lattice is used as an initial guess for the algorithm. To fit data, one assumes the longitudinal position of the beam’s magnets are fixed and known. BPM errors and the precise strength of each magnet are then treated as variables and a numerical algorithm is iterated to fit the clean BPM data with those variables. The machine’s ideal lattice is used as an initial guess for the algorithm. To fit data, one assumes the longitudinal position of the beam’s magnets are fixed and known. BPM errors and the precise strength of each magnet are then treated as variables and a numerical algorithm is iterated to fit the clean BPM data with those variables. The machine’s ideal lattice is used as an initial guess for the algorithm. To fit data, one assumes the longitudinal position of the beam’s magnets are fixed and known. BPM errors and the precise strength of each magnet are then treated as variables and a numerical algorithm is iterated to fit the clean BPM data with those variables. The machine’s ideal lattice is used as an initial guess for the algorithm.

The data on BPM gains and tilt from several sets of data can be analyzed for statistical patterns. Such an analysis was done on the BPM error data for April 29th, 2003 (pre-processed by MIA), and the results show that some BPMs have consistent systemic errors. These errors, along with data about noisy and bad BPMs from the FFT and SVD procedures, are being presented to the PEP-II BPM Working Group with a reccomendation for recalibration.

While the scope of this project is limited to pedagogical exploration of linear optics due to time constraints, I strongly hope to continue this project with the purpose of using MIA to find an analytical method to determine the Courant parameter and betatron oscillations in the eigenplanes of a coupled system.

References

Woods, H. R. (2003).
PEP-II Attains Record Beam Luminosity. Interaction Point.
May 16, 2003.

2003-07-24T15:04:00.000-07:00

I isolated my ringing-in-the-middle BPMS along with my noisy BPM. And I'm having a lot more fun plotting stuff on matlab! Things are fun again.

I would still like to study the frequency shizzle, however. I wonder if I can find someone who can help me with that. (How to interpret frequency content)

Maybe i'll look at that old printout I had.
Speaking of which--reminder: get old printouts and materials.

2003-07-24T14:18:00.000-07:00

Went back and finished an old project: custom 4x4 ideal lattice creation. It might come in handy some time in the future. I could always customize a design to form an "idealized" linear PEP-II with no energy consideration. Of course... then it would be the same as any other machine of similar idealization.

Back to those fourier transforms.

2003-07-24T11:47:00.000-07:00

I'm frustrated by oscillating between believing that I have nothing to write about, and having everything in the world to write about. I'm either writing an article for the Stanford Daily or an encyclopedia.

2003-07-24T11:06:00.000-07:00

Here's something cool that I forgot to mention from yesterday:

PEPII uses RF frequency microwaves to regulate the beam and keep it "bunched up." It is conventional to think of the waveguides as slowing down the electrons ahead of the bunch and speeding up those behind it.

While this is the net effect, what actually happens is that the electrons ahead of the bunch are sped up, causing relativistic effects and making the particle more massive. Massive particles have more intertia and do not turn as tightly, and thus are relegated to a slightly larger (radius) orbit. The size of this larger orbit overcompensates for the higher speed, and the particle falls behind to join the rest of the bunch.

Pretty neat, huh?

My questions about what i'm doing.

MIA DATA
(1) BPM gains and tilt. Easy to compute, simple statistics.
(2) BPM data
(a) Identify bad, dead BPMS
(b) Extract clean signal
(c) Find eigenplane

Questions:
How does this all lead to higher luminosity?
What is my overreaching goal?
So even if i can make a linear optics representation--which is only an academic contribution--what is really done to modify the machine?
What is nonlinear optics?

2003-07-23T16:40:00.000-07:00

Assignments:
* Use SVD (After eliminating closed orbit) and compare results to FFT

NOTE: don't worry too much about getting rid of the funky spikes on the end, just use the spikes that are correct

* try using FFT and SVD together

Some values:
range of betatron oscillations: .5-.7
range of synchrotron oscillation: .03

Food for thought: figure out what the freakin' units are on an FFT'd function. They should be something like __*2pi = angular frequency.

Long term project:
Can I extract eigenplane 1 and 2 data from bpm data? (This is publishable work)

2003-07-23T14:59:00.000-07:00

What do I do with the two dead BPMs that I found?

2003-07-23T14:51:00.000-07:00

I can do it--I can do it--I can do it--I can do it--I can do it--I can do it!!

I have to filter by coinciding peaks first!!
Then I can make a vector of coinciding peak and then use the SAME "chop of at a certain discrimination fraction" technique that I've been doing all this time.

2003-07-23T10:16:00.000-07:00

Yesterday: Talk by Andrei Linde. He's an *excellent* speaker and absolutely hilarious.

Today, goals:
1. program to take noise out of bpm data
2. program to convert clean frequency content into clean signal
3. analyze clean signal
4. SVD enhance
5. analyze enhanced signal

2003-07-22T15:40:00.000-07:00

Rahul just taught be a little about convolution (and he recommended another book for me to learn FFT from). Very cool stuff!

2003-07-22T14:14:00.000-07:00

Observations from playing with april 29 bpm data:
Funny xbuferX results for BPM 1 (xcomTags(1)) and 318 (xcomTags(196)).

Why are they at such a different range?

2003-07-22T11:34:00.000-07:00

Finished the main part of my BPM error program. I think I'll go back to playing with FFT and the BPM data. I wonder what else there is for me out there?

2003-07-21T15:43:00.000-07:00

Just got the April 29 Data
My goals:

1. FFT, isolate closed orbit (see notes on board about what a closed orbit is--it has to do with static BPM shifts w/rt the beam)

2. SVD, try plotting U and V'

3. See what I can do with the orbits. Can I make an analysis of the BPMs independent of the model?

4. play with FFT.

Here's something to think about: M12 = A2*R12*A1^-1

2003-07-21T13:54:00.000-07:00

Need to develop a test program for my runAnalysis program... dataArray matrix variabl especially. I'm thinking either a full out '==' comparison or using random integers.

2003-07-21T13:09:00.000-07:00

personal:

Talked to Dr. Yan briefly about long-term opportunities working under him. I'm curious where this might all lead up to. I should also remember to talk to Scott Thomas and maybe Pat Burchat/John Fox/and company. Yue recommended talking to Dr. Peskin, her graduate advisor.

I should go ask Erin Smith about accessing previous students' work.

Flip loves physics.

2003-07-21T12:24:00.000-07:00

More books from the library:

Introduction to Lie Algebras and Representation Theory, Humphreys
Application sof Lie Groups to Differential Equations, Oliver
Introduction to Statistical Analysis, Dixon
The Art of Electronics, Horowitz and Hill

Borrowed from John Fox:
Student Manual for the ARt of Electronics
Electronic Circuites, Tietze
Principles of Electronic Instrumentation, Diefenderfer

2003-07-21T09:32:00.000-07:00

A few notes:
* Working on BPM gain/coupling errors--I should develop a program to analyze it... time to dig up some stat
* On the back burner: should finish up lattice program--results should be interesting
--goal: input any lattice (by hand or by text), create all appropriate matrices, output graphs of betatron function, output all releveant frequencies using fft. This will be a pretty cool program!

Personal notes:
* Should probably purchase Horowitz & Hill + Student Manual
* I think I'll also look into the "teach yourself electronis" book
* I need to return/renew my library books (note: Rahul still has "6 not-so-easy pieces")

2003-07-18T11:20:00.000-07:00

Developed a testfft script that clarifies a lot of things. I will soon post it on my website.

echo on
% Test FFT Script
% July 18, 2003 by Philip Tanedo
% Purpose: demonstrate the basic use of matlab's fft function
% Adapted from: "Spectral Analysis using the FFT" by Brett Ninness
% ... http://www.ee.newcastle.edu.au/brett/elec2400/matlab4.pdf

ts=0:.1:50;
% Time data (independent variable)
% The only really important piece of data here is that dt=.1
% note: we usually start with x, not ts

x=sin(ts)+cos(3*ts);
% Construct sampling vector
% Usually, this is our starting point (and we know what dt is)

fs=10;
% Frequency = 1/t (frequency of sampling)

t=(0:1:length(x)-1)/fs;
% Construct time vector
% This is the same thing as ts, but we construct it now
% because we don't usually start with ts

figure
plot(t,x)
xlabel('time, .1 second intervals')
ylabel('x(t)')
title('Plot of x(t)')

N=length(x); % For convenience

ws = 2*pi/N;
wnorm=-pi:ws:pi;
wnorm=wnorm(1:N); % Because length(wnorm) != length(X)

w=wnorm*fs; % scaling

X=fft(x);
Xs=fftshift(X);

figure
plot(w, abs(Xs))
xlabel('frequency, rad/s')
ylabel('magnitude')
title('Plot of abs(fftshift(fft(x(t))))')

figure
plot(w, real(Xs))
xlabel('frequency, rad/s')
ylabel('magnitude')
title('Plot of real(fftshift(fft(x(t))))')

figure
plot(w, imag(Xs))
xlabel('frequency, rad/s')
ylabel('magnitude')
title('Plot of imag(fftshift(fft(x(t))))')

echo off

2003-07-18T10:25:00.000-07:00

Yesterday I talked to John Fox about Applied Physics at Stanford and he pointed me towards a Caltech freshman lab class implementing the ZAP! take-home labs.

Also, Persis Drell gave the best talk that our group has had. It wasn't about physics (surprisingly)--it was about her life and what it's like to *become* a physicist. It was very entertaining and very up front.

I'm continuing to work on FFT. It turns out you need to use fftshift to make these stupid plots look right. I'm still trying to figure out how to scale them properly.

2003-07-17T13:20:00.000-07:00

John Fox is always busy.

FFT, IFFT.

2003-07-16T12:40:00.000-07:00

Cleared up a lot of questions.

My assignment:
(1) Understand
(2) FFT/SVD on data

Questions
1. What do the sine and cosine components of FFT *mean*? Maybe Weaver can help me. I just want to make sure that I have the concept down pat. What about the imaginary/real parts?

2. Make sure I understand the paper I was reading.

Comments
1. Q is essentially guessed using bpm data. Q is in eigenspace and is guessed using real space data. (Everything else is in real space)
2. The (for example) 173 bpms are used *TOGETHER* to form a cohesive picture of what's going on... also to figure out other elements of the R matrices. Consider: Rab, where only Rab12,14,32,34 are known. Ditto with Rbc. We can also construct Rac knowing Rac=RbcRab. We now have information about the other terms in Rac. Magic!

2003-07-15T15:56:00.000-07:00

M'SM=S
Zf=MZi

M=Zf*inv(Zi)
Zf'SZf=Zi'SZi=Q
inv(Zf')*Q*inv(Zf)=inv(Zi')*Q*inv(Zi)=S

inv(Zi)=inv(Q)*Zi'*S
Which leads (naturally, though I didn't grind the numbers) to equations 1-4 in "Application of Model-Independent Analysis to PEP-II Rings", SLAC-PUB 8888.

I'm still not sure what space M is in. I believe we're talking about eigenspace, but then Z would be terribly boring. Also, then M14 and M32 are zero.

I'm confused by the latter part of section 3 of the paper where Dr. Yan discusses BPM gains and x-coupling multipliers.

2003-07-15T11:20:00.000-07:00

Flurry of new ideas that I need to get a handle on:
* Green's functions
* "Local transfer map"
* More ODE stuff

* LaTeX

2003-07-14T16:45:00.000-07:00

My new project: http://accelconf.web.cern.ch/AccelConf/p01/PAPERS/RPPH134.PDF

2003-07-14T15:42:00.000-07:00

IT WORKS!

2003-07-14T14:10:00.000-07:00

I found the factor of square root of 2. Congrats, Flip. Congrats.

2003-07-14T12:14:00.000-07:00

Once again reluctance is the bane of my existence.
I should be more like Max Planck. He derived a formula for blackbody radiation that worked, and THEN set out to figure out what it meant. I was too reluctant to use my /\ matrix, and it turns out it was correct.

Well, it was correct in theory. I applied it incorrectly and got the wrong result. I can pat myself on the back for having good intuition. Too bad I didn't go all the way with it.

The lesson: We have M. We want R. We don't care what A is. The A matrix is not sacrosanct. If I have to rotate, twist, transform it in any arbitrary way to make it symplectic, then so be it.

I should finish my program within the hour.

2003-07-14T11:06:00.000-07:00

The rest of my program is written. I still need to make a symplectic Att.

2003-07-14T09:51:00.000-07:00

Nope. /\ only works when s=0. A very special and unrealistic case.

2003-07-14T09:37:00.000-07:00

A weekend has not solved my problem, but it has calmed me down a little bit.
I think there may be some hidden constraint that I have overlooked.

Thus far I know that M is symplectic. I know that fastDecompose (specifically the decomposeM function) works the way I want it to and that At is of the correct format.

But here's the thing: why should At be symplectic with respect to a scalar gauge? I don't think it does. Even when s=0, the X and Y invariant subspaces are transformed differently by M. (From the beginning: each quad is focus/defocus) As At is a matrix of real and imaginary parts of the complex eigenvectors ((note: i still want to figure out the constraints on s)) and as the two subspaces of M are different, then it is clear that even when s=0, At wouldn't have a single gauge scalar.

I will ask dr. Yan about this later. For now I will use /\= [Sc1*I 0; 0 Sc2*I] where each element is a 2x2.

2003-07-11T17:29:00.000-07:00

This is definitely very frustrating. I'm heading out 10 minutes early so I can get a head start on jogging on the Dish. Maybe the cows would know something or another about symplectic matrices.

2003-07-11T17:19:00.000-07:00

Just got started with my new approach ... and another hitch.
I'm not getting a consistent scalar gauge.

Troubleshooting:
* M matrix IS symplectic. So at least everything up to there is correct. (Lattice design)
* Decomposition function is correct (i.e. it does what I want it to do. Whether or not what I want it to do is correct or incorrect is another question)

But my A-tilde matrix (At) cannot be made symplectic by a scalar gauge. Why?
There are only two lines of code between the M matrix and At matrix.

[V,E] = eig(M) % makes V the matrix of eigenvectosr
At = [real(V(:,1)) imag(V(:,1)) real(V(:,3)) imag(V(:,3)) % Assigns At as I wanted

I perform At'*S*At and I get a diagonal block matrix of 2x2s. But the upper left block does not agree with the lower right block, and hence two different gauges exist for each subspace--which is NOT mathematically sound: there can only be one gauge.

I wonder if my A-tilde is defined correctly.

2003-07-11T16:10:00.000-07:00

Wait. Maybe it doesn't have to be symmetric.
The quadrupoles aren't symmetric, so why should a particle's one-turn map (M) be symmetric with respect to the invariant eigenplanes (s=0)? I'll have to try plugging in -k for k in my lattice and then checking if the resulting M is the same thing as M....

OH WAIT. A visit from Dr. Yan:
I was doing everything in the wrong order!!

(1) The answer to the above question is yes.
(2) Here's my modus operandi.

(a) [R,At,M]=fastDecompose(k,L,s) (this step is correct)
(b) transform At -> Att = Sc*At such that Att is SYMPLECTIC (AHA!!!)

note: Dr. Yan TOLD me this before and I didn't understand its significance until now. My problem was that I put the scalar and rotational gauges together. No! The scalar gauge does one thing (makes Att symplectic) and the rotational gauge does something completely different (rotates into "proper form").

(c) Find Rt (My decomposeAx and decomposeAy functions were completely incorrect).
(d) find (CA)
(e) perform my end-analysis to find C and alphas, betas.

2003-07-11T14:22:00.000-07:00

Something *is* indeed wrong. I can't believe that I didn't check this earlier, but I do *NOT* end up with a symmetric matrix when I put in skew = 0!!

Where is the error?

2003-07-11T14:07:00.000-07:00

I noticed that all my Rt matrices are identity. I will modify my program to allow for any arbitrary lattice design. Maybe the problem lies in the skew matrix that has no counter skew?

To do:
* backup semi-complete 4x4coupled folder
* develop new createLattice program

2003-07-11T13:59:00.000-07:00

I already found one error in the decomposeAy function. (It was written decomposeAx)

Still the same errors, though.
The ScC matrix is totally incorrect. The two diagonal 2x2 submatrices do not have zero off-diagonals. Meanwhile, the two off-diagonal 2x2s don't have matching Ws.

2003-07-11T13:57:00.000-07:00

This morning we went to go visit Applied Materials in Santa Clara.

I finished scripting my program so that all I have to do is run a single command line prompt and enter a few values at the beginning. The problem is, I'm still getting incorrect scalar-times-C matrices. I'm wondering if there is *another* gauge rotation involved, but I don't think that would fit into the framework of M=At*R*At^-1

I typed up all of my functions and scripts. I include them below for no apparent reason.

--- (i wonder if the below will turn into a line)

% initialize k,L,s

clear all;
global k;
global L;
global s;
initializekLs;

function initializekLs
% function initializekLs
% preconditions: k,L,s global variables; no other variables
% postconditions: those variables filled with ik,iL,is
% use with initialize.m script

global k;
global L;
global s;

k=input('Input k: ');
L=input('Input L: ');
s=input('Input s: ');

% algorithm to partially solve 4x4
%
% INPUT: k,L,s
% fastDecompose: gives R,At,M
% decomposeAx: gives Alx,Bex,Thx
% decomposeAy: gives Aly,Bey,Thy
% createA: gives A
% createRt: gives Rt
% createScC: gives ScC
% *one more function to separate Sc, Phi, and W
% *--not yet constructed
%
% Preconditions: k,L,s already exist
% Postconditions: all variables filled

[R,At,M]=fastDecompose(k,L,s);
[Alx,Bex,Thx]=decomposeAx(R,At,M);
[Aly,Bey,Thy]=decomposeAy(R,At,M);
[A]=createA(Alx, Bex, Aly, Bey);
[Rt]=createRt(Thx, Thy);
[ScC]=createScC(At,Rt,A);

%TESTING: view ScC
ScC

function [R,At,M] = fastDecompose(k,L,s)
% function [R,At,M] = fastDecompose(k,L,s)
% [R,A]=decomposeM(createM(createLattice(k,L,s)));
% THIS FUNCTION DOES NOT WORK for LARGE s

[R,At]=decomposeM(createM(createLattice(k,L,s)));
M = createM(createLattice(k,L,s));
-------------------------------------------------------
function [lattice]=createLattice(k,L,s)
% function [lattice]=createLattice(k,L,s)
% creates 4x4 coupled lattice
% k = focusing
% L = drift
% s = coupling
% coupling occurs after o1
% map: o1 - s - o2 - o3 - o2
% o1 = x focusing, y defocusing
% o2 = drift
% o3 = x defocusing, y focusing
% s = coupling
% lattice = [o2, o3, o2, s, o1]
% note convention: M = L(1)*L(2)*L(3)*L(4)*L(5)

o1 = [1 0 0 0; -k 1 0 0; 0 0 1 0; 0 0 k 1];
o2 = [1 L 0 0; 0 1 0 0; 0 0 1 L; 0 0 0 1];
o3 = [1 0 0 0; k 1 0 0; 0 0 1 0; 0 0 -k 1];
S = [1 0 0 0; 0 1 s 0; 0 0 1 0; s 0 0 1];

lattice(:,:,1)=o2;
lattice(:,:,2)=o3;
lattice(:,:,3)=o2;
lattice(:,:,4)=S;
lattice(:,:,5)=o1;

function [M]=createM(lattice)
% function [M]=createSimpleM(lattice)
% creates M as a composition of 4x4 matrices in the lattice
M=lattice(:,:,1)*lattice(:,:,2)*lattice(:,:,3)*lattice(:,:,4)*lattice(:,:,5);

function [R,At] = decomposeM(M)
% function [R,At] = decomposeM(M)
% M = AtRA^T-1
% At = A-tilde (At = CA)
% C = [ Icos Wsin] A= [ Ax 0 ]
% [-Wsin Icos] [ 0 Ay ]
% I, W, Ax, Ay are 2x2
% Ax defined conventionally as [ sqB 0 ]
% [ a/sqB 1/sqB ]
[V,E]=eig(M);
At=[real(V(:,1)) imag(V(:,1)) real(V(:,3)) imag(V(:,3))];
R=inv(At)*M*At;

function [Alx,Bex,Thx]=decomposeAx(R,At,M)
% function [Alx,Bex,Thx]=decomposeAx(R,At,M)
% decomposes Atx into Ax, where Ax is the nice decoupled form
% by giving values for alpha-x (Alx) and beta-x (Bex)
% also gives value for the x-rotation Thx
%
% Does NOT give a value for phi or Sc
% But I know Sc = 1/det(At)*cos(phi)

% some convention
% Atx = [ m n ] = top left 2x2 of At
% [ o p ]
m = At(1,1);
n = At(1,2);
o = At(2,1);
p = At(2,2);

Atx=[m,n;o,p];

% My system of nonlinear equations
% (0) Sc^2*det(Atx) = cosPhi
% (1) Sc(mcosThx-nsinThx) = cosPhi*sB
% (2) Sc(ocosThx-psinThx) = -Al/sB *cosPhi
% (3) Sc(msinThx+nsinThx) = 0
% (4) Sc(osinThx+pcosThx) = 1/sB *cosPhi

% from (3)
Thx=atan(-n/m);

% more convention
a1=(m*cos(Thx))-(n*sin(Thx));
a2=(o*cos(Thx))-(p*sin(Thx));
a3=(m*sin(Thx))+(n*cos(Thx));
a4=(o*sin(Thx))+(p*cos(Thx));

Bex=((1/det(Atx)*a1)^2);

% from (1) & (2)
Alx=-(sqrt(Bex)*a2)/(det(Atx));

% from (0)
% Sc^2 =(1/det(Atx))*(cos(Phi))^2;

function [Aly,Bey,Thy]=decomposeAy(R,At,M)
% function [Aly,Bey,Thy]=decomposeAx(R,At,M)
% decomposes Aty into Ay, where Ay is the nice decoupled form
% by giving values for alpha-y (Aly) and beta-y (Bey)
% also gives value for the y-rotation Thy
%
% Does NOT give a value for phi or Sc
% But I know Sc = 1/det(At)*cos(phi)

% some convention
% Aty = [ m n ] = top left 2x2 of At
% [ o p ]
m = At(3,3);
n = At(3,4);
o = At(4,3);
p = At(4,4);

Aty=[m,n;o,p];

% My system of nonlinear equations
% (0) Sc^2*det(Aty) = cosPhi
% (1) Sc(mcosThy-nsinThy) = cosPhi*sB
% (2) Sc(ocosThy-psinThy) = -Al/sB *cosPhi
% (3) Sc(msinThy+nsinThy) = 0
% (4) Sc(osinThy+pcosThy) = 1/sB *cosPhi

% from (3)
Thy=atan(-n/m);

% more convention
a1=(m*cos(Thy))-(n*sin(Thy));
a2=(o*cos(Thy))-(p*sin(Thy));
a3=(m*sin(Thy))+(n*cos(Thy));
a4=(o*sin(Thy))+(p*cos(Thy));

Bey=((1/det(Aty)*a1)^2);

% from (1) & (2)
Aly=-(sqrt(Bey)*a2)/(det(Aty));

% from (0)
% Sc^2 =(1/det(Aty))*(cos(Phi))^2;

function [A]=createA(Alx, Bex, Aly, Bey)
% function [A]=createA(Alx, Bex, Aly, Bey)
% creates 4x4 A matrix in proper form

Ax=[sqrt(Bex) 0;-Alx/(sqrt(Bex)) 1/(sqrt(Bex))];
Ay=[sqrt(Bey) 0;-Aly/(sqrt(Bey)) 1/(sqrt(Bey))];
A=[Ax, zeros(2); zeros(2), Ay];

function [Rt]=createRt(Thx, Thy)
% function [Rt]=createRt(Thx, Thy)
% creates a 4x4 gauge rotation matrix
% NOTE: this actually produces Rt^-1
% st ScAt(Rt^-1)=CA

Rx=[cos(Thx) sin(Thx);-sin(Thx) cos(Thx)];
Ry=[cos(Thy) sin(Thy);-sin(Thy) cos(Thy)];

Rt=[Rx zeros(2); zeros(2) Ry];

function [ScC]=createC(At,Rt,A)
% function [ScC]=createC(At,Rt,A)
% creates scalar*C by: ScC = At*Rt*A
% note: Rt is really Rt^-1
% (in all my scripts we never use the "proper" Rt
% so for efficiency, I write Rt^-1 as Rt.

ScC=At*Rt*A;

2003-07-10T15:39:00.000-07:00

Two new functions: createA and createRt. I tested these by producing At*Rt*A^-1, which is supposed to produce a scalar multiple of C, but the results came out blatantly incorrect.

I'm not sure where to start troubleshooting now. Perhaps its time to call it a day. We have a lecture at 4:30 by Paul Kunz, and I'd like to purchase my Quantum Mechanics book at the bookstore before then.

2003-07-10T15:03:00.000-07:00

Hm. Plugged in scalar gauge (I've been careless with notation... I refer to it as Sc here, and write it as either a lambda or lambda-tilde in my notes) and that made my system redundant--i.e. I have a relation between lambda and cos(phi), but right now that relation is arbitrary. Perhaps solving for the C matrix will clear this up.

Interesting developments.

Status: I can now (I think) extract alpha and beta from R,At,M. I also have an extra equation in my system that I can use to self-check. I can no longer find values of phi or lambda... which means I can no longer check for correctness between the invariant X and Y subspaces. But oh well--that's what I "intuitively" expected, in the first place.

2003-07-10T14:22:00.000-07:00

Aside: I would like to purchase some introductory electronics kits so I can learn about circuit analysis and do all that fun stuff. I'm willing to invest and entire paycheck towards the cause of learning electronics.

2003-07-10T14:13:00.000-07:00

AHA! (1) no, Rt must be coupled. ((Otherwise what's the point of calling C the decoupling matrix?)) (2) I found an error that might fix things: I forgot to include the scale gauge in my calculations!

Lesson: If there is no immediately apparent error in the program, write out the equation by hand enough times until you realize you were neglecting one of the elements.

2003-07-10T13:42:00.000-07:00

Potential solution: I think I oversimplified when I assumed that Rt^-1=[Rx,0;,0,Ry]. R is a general gauge rotation, perhaps it is coupled?

2003-07-10T13:19:00.000-07:00

A problem: I ran my lattice thus far with (.3,.4,.03) and got values of phi that were different for the decomposeAx and decomposeAy functions. Something is fishy here.

2003-07-10T13:13:00.000-07:00

A complete solution is not far away.

2003-07-10T10:29:00.000-07:00

Oh! Now it's easy! I just need to figure out MatLab's rules for inverse trig functions, but I can now find all my unknowns for the x phase space. It should be even easier in the Y phase space.

Is the gauge scalar the same for the X and Y subspaces? It must be.

2003-07-10T10:24:00.000-07:00

Retroactive note for yesterday (Wed, Jul 9): Last night we also had the first class of the Travis School of Quantum Mechanics. We discussed the 1-dimensional wave equation. Something completely off-topic that I learned: I can derive trigonometric identities *VERY* easily by using Euler's equation. Fascinating!

Currently: I'm working on the AtRt^-1=CA equation that I was working on yesterday. I've reached the "ugly" step which involves a nonlinear system of equations requiring some "clever algebra."

Some discoveries: let the top left 2x2 of At (At is the eigenvector-derived matrix) be At'' (A-tilde 1,1).
Sc*det (At''Rt^-1) = cos^2(phi)det(Ax) = 0, where Ax is the nice form of the decoupled A matrix for the x phase space.
Thus, as the det of a rotation is 1 and the det of the composition of matrices is the product of their individual determinants:
Sc*det(At")=cos^2(phi)

I now have a way to find the decoupling angle given the scalar gauge and vice versa.

When I solve the equations, I end up with a nonlinear system (sines, cosines, square roots... oh my). I will continue to work and see if I can get a nice, analytic solution.

2003-07-10T08:55:00.000-07:00

Books Checked Out From SLAC Library:
Applications of Discrete and Continuous Fourier Analysis, Weaver
Fourier Analysis and Generealized Functions, Lighthill
Introduction to the theory of Fourier's series and integrals, Carslaw
Problems and Solutions in E&M
Quantum Physics: An Introduction
Differential Forms, Schrieber
6 Not-So-Easy Pieces
(Additionally, there was some other book on logic and ideas... but I returned that the day I checked it out)

2003-07-09T17:41:00.000-07:00

FFT
X(k) = sum x(n)*exp(-i*2*pi*(k-1)*(n-1)/N), 1 <= k <= N. ((sum from n=1 to N))
(from matlab)

Dr. Yan showed me an example of how to use FFT on some data that we made up on the spot. I still need to learn the theory before I can really grasp the mechanics of the fourier transform.

Insight on my current goal: AtRt=CA
I want to extract alpha-x, beta-x, and alpha-y, beta-y and ultimately find C.

At the A matrix extracted from the eigenvector matrix of M.
Rt is the yet-unknown gauge rotation. (form: Rt = [Rx,0;0,Ry], each a 2x2)
C is the yet-unknown decoupling matrix (known form, 5 unknowns)
A is the A matrix in proepr form ([Ax,0;0,Ay]).

I came up with this on my own from straightforward analysis... but I foolishly gave up on the approach because I thought there were too many unknowns and not enough constraints. (As it is, that may yet be the case.) However, Dr. Yan suggested the very same approach (not knowing I already had it outlined in my lab notebook), and pointed out some subtleties (i.e. the existence of Rt--there is always a gauge transform!) that make my life a lot easier.

The idea is that the diagonal 2x2 elements of C are identity matrices multiplied by a cosine. This means the diagonal 2x2 elements of (CA) are just the diagonal 2x2s of A multiplied by the cosine. Thus, I can treat this as nearly the exact same thing as the 2x2 case! There are, of course, some complications: I have a gauge rotation Rt (which is also decoupled) and the angle of the cosine needs to be determined. All in all, that is 3 extra constraints. I now have, however, two decoupled systems. I believe that with some quality time spent in the library, I should be able to extract alpha and beta for x and y from this information.

Once alpha and beta are extracted, things begin to unravel very nicely. I will then know At, Rt, and A. C = AtRtA^-1. From there I can determine cosine ... hey, I think that's what Dr. Yan was trying to explain to me: In the above case I don't think I need to figure out what the cosine is, only that it is the same for the x and y systems! Anyway, from there I can determine cosine and the four elements of W by clever algebra.

Never trust someone who uses the phrase "clever algebra."

This is very exciting. Maybe my computer will have fixed itself and today will be a near-perfect day.

2003-07-09T14:08:00.000-07:00

Can I extract alpha and beta from R and M in the 4x4 coupled case?

2003-07-09T13:36:00.000-07:00

4by4coupled program works fine, actually.
It turns out I was just being irresponsible with my value for s.
s << k but, by what factor? this is something I am very curious about...

In the 2x2 case, I was able to use an important fact:
M^TSM=S <==> det(M)=1
((Thus far, I have found no simple "check" for symplecticity in the 4x4 case.))
Using this equivalence, I was able to explicitly find the bound on k and L in a simple lattice (similarly, I could find the bound for any other 2x2 lattice) using det(M)=1, instead of the nebulous symplecticity definition.

I wonder what I can extract from 4x4 symplecticity that would be useful in providing constraints on s for a given lattice. All I know right now is that s needs to be small.

Also notable: it turns out the matrix that I derived yesterday has a special name... some kind of Courant matrix. Dr. Yan was trying to show it to me and I said "oh yes! I have that somewhere in my notes." Of course, it only took him 2 minutes, whereas I spent a good 30 minutes deriving it on my own.

2003-07-09T12:26:00.000-07:00

A pedagogical breakthrough! The eigenplane is the coordinate system relative to the beam (i.e. the betatron oscillation in the eigenplane is decoupled). The eigenplane wobbles with the skews.

z st z(i+1) = M(zi)

then Z = matrix of 4 such z's (to be specified): Z=[z1 z2 z3 z4]

Z(i+1)=MZ(i)
M = Z(i+1)Z^-1(i)

What are zi's? Correspond to four degrees of freedom for orbit. The axes (I think... or something along those lines) of the eigenplanes.

2003-07-09T11:28:00.000-07:00

Posted new web blog to accompany my lab book.

I'm still trying to figure out how to extract C, the decoupling matrix given M and R (and A for Md, where Md is the non-skewed matrix). I thought I would just use MC=CMd, but these calculations ended up becoming tedious and the results were ultimately trivial--i.e. I ended up with more statements about constraints on Md and M than I did about the unknowns. (Although maybe I just didn't look hard enough)

There must be an easier way to extract C.

2003-07-09T10:29:00.000-07:00

A recap of what I've done thus far:

Monday, June 23
-Lecture: Introduction to SLAC, Helen Quinn
Basic "what's a meson" introduction to particle physics

-SLAC Tour
I wish I had brought my camera

-Lunch with Mentors
I met Dr. Yiton Yan today. I also met Rahul, the EE grad student (for Martin Lee) sharing my office with me.

-Worksite with Mentors
Read Dr. Yan's lecture notes on Lie Algebraic methods for beam accelerator physics.

Tuesday, June 24
-Safety and GERT training
I felt pretty bad about getting drowsy during this... but I learned a few cool things: the "powers that be" (FDA or some similar agency) say that an annual dose of 5000 mrems of radiation is "safe." Americans receive an annual dose of 360 mrems. As a GERT employee, I am allowed 100 mrems of SLAC-related radiation. I'm not supposed to go and spend it all in one place.

-Got my badge and dosimeter
It was like a prize for staying awake through GERT training.

-Visited the control room
Dr. Yan was very hopeful that they would be able to bump up the luminosity to a new record (MIA set the record at 6.2 cm^-2 s^-2 earlier this year), but there were complications with the machine. These were later cleared up, but futher errors (human and ideological) aborted the run. Here's what was interesting and edifying, though: usually the accelerator engineers/physicists who dial in the magnet modifications make drastic changes in percentages... i.e. instead of shifting all the magnets all the way, all at once; they shift the magnets each by a fraction of the total change and then make sure the beam is still stable. An unstable beam would be a catastrophic error. Anyway, Dr. Yan noticed that the current set up of the machine was very noisy--the betatron oscillations were kind of wild... but even then the luminosity was very high. He was hoping that by using MIA to tame the betatron oscillations, he would be able to bump the luminosity even higher. The noisy beta, however, led MIA to suggest a drastic solution. When the engineers were dialing in the solution, the beam destabilized, and the run was aborted. Why? MIA's solution was nonlinear, and they were dialing in the solutions linearly (i.e. 10% of the solution across the board... 20% of the solution, etc.).

Dr. Yan mentioned this possibility after the run was aborted. In the meanwhile, I learned some basic MatLab commands for MIA analysis.

Wednesday, June 25
Dr. Yan was right. He analyzed the step-by-step dialing of solutions, and he found that even on the virtual machine the beam would destabilize using that method. He used MIA to develop a nonlinear, incremented solution... but unfortunately, he was unable to ever implement this because the beam shuts down after this week and the machine will be changed.

I continued to work on M=ARA^-1, looking for an analytical solution showing that this was implied from S=M^tSM, the symplectic condition. No real breakthroughs, but I did some nice calculations by hand. Rahul taught me a little bit about SVD, and that was very fascinating. I'll have to review it later.

Today I also met Rahul's mentor, Martin Lee. He's a very fun guy--he laughs a lot and makes the room a lot livelier when he's around.

Thursday, June 26
3pm: Dr. Yan presents his findings.

I continue working on my symplectic<==>normalized rotation proof in the 2x2 case. Dr. Yan hinted at using eigenvectors. At the time this was very nonintuitive... but reading Bretscher made this much clearer.

I also took a break today to go to the library. I check out some books and thought about Maxwell's Demon and entropy (the weather is conducive to that).

Showed Dr. Yan my first draft of my project description.

Friday, June 27
My goal was to get this M=ARA^-1 thing DOWN by today. And I think to a reasonable extent, I did. I was also able to get my computer account set up.

Took notes on Bretscher's Linear Algebra text on SVD. The process makes much more sense. It's an approximation method... kind of like a Taylor series.

Finalized Project Description.

What I know thus far:
-for a 2x2: symplectic <==> det =1. unfortunately this does not hold for the 4x4 case.
-for a symplectic M that can be reduced to M=ARA^-1, 0<(a+d)^2<4. This didn't make too much sense to me initially, but Dr. Yan commended me on the discovery and said that if this condition did not hold, then A would not be invertible and M would be an unstable one-turn-map. It makes sense!
-So here's how it works: we get M by analyzing the lattice elements (quad matrices, K, and drift matrices, L)... then we decompose into ARA^-1. I later (one week later) found out that we decompose this FURTHEr into A = ScAtRt, where Sc is a scalar, At is an "A" matrix in "proper form" (with alpha and beta), and Rt is a gauge rotation. This is fine because we en dup with: M=ScAtRtR(1/Sc)Rt^-1At^-1 = AtRtRRt^-1At^-1... but the three middle matrices are rotations, so this is equivalent to M = AtRAt^-1. Thus, At is a valid A.
-For a 4x4, this is a little more complicated and there is a decoupling matrix, C. M=CARA^-1C^-1.

Monday, June 30
First day of lecture week. The beam is down and most of the researchers aren't going to work this week.

-Lecture: How Accelerators Work, John Fox
Decent intro... I wasn't compelled to take very detailed notes (the slides went by very quickly and his topics were a little haphazard). Once again Microwave-dude from QuarkNet makes a fuss about... microwaves.

-Accelerator Tour, John Fox
Once again I forget my camera. I'm amused by all of the photos of people pretending to press "emergency off" buttons and such.

-Lecture: Particle Detectors, Colin Jessop
This was actually a decent lecture--but I learned a lot about how sleepy I was. (And the rest of the group... Helen Quinn included) I later ended up looking up Dr. Jessop's website to learn a little more about the things that I missed by being drowsy... to little avail. I found most of the material in Helen Quinn's book (The Charm of Strange Quarks), however, and I no longer feel like I missed out.

Tuesday, July 1

-Lecture: Quantum World and Particles, Helen Quinn
-Lecture: Hidden Scatterer/Feyngames, Helen Quinn
I now have a good understanding of how to use Feynman diagrams and my Particle Data Group book. A very happy morning!

-Exploration: Cosmic Rays
Learned about cosmic ray scintillator detectors.

Wednesday, July 2
-Lecture: Dark Matter and GLAST, Larry Wai and Eduardo do Couto e Silva. Dr. Wai's lecture was a little dense (well, both were)... I need to read more about Majorana and Dirac particles. GLAST is something Sandi might know a little about.

Thursday, July 3
-Lecture: BABAR-What For? Pat Burchat
A *GREAT* lecture from my future P134 professor. See notes for good stuff about CP violation.

-Eduardo poses a riddle: Why do pions decay preferentially into muons, not electrons? (I found the answer about a week later: it has to do with conservation weak isospin and the fact that leptons have mass. Counter-intutively, the particles' masses are what allows them to violate weak isospin. It doesn't decay into taus often because taus are too massive.)

Monday, July 7
Back to work... put all my 2x2 discoveries into MatLab.
I did more analysis on my simple lattice. Found a few things.
-calculated an explicit M for given k,L
-conditions on k,L s.t. M is symplectic and M=ARA^-1: k^2L^2<2. (This makes creating 2x2s a lot easier using my MatLab program... now I know why I kept getting spurious results for "ugly" k-L combinations.)
-Tried to find a non-zero closed one-turn orbit... does not exist. (Proved this)

Instructions: start working on a 4x4 with skew quads
-found 6 constraints on a symplectic 4x4... not immediately useful to me. (They are related to DoF of the hamiltonian)
-found that KS = SK.

Tuesday, July 8

How is C related to S?
I have found that the diagonal 2x2 submatrices of any M whose lattice has a skew matrix are preserved...i.e. are equivalent to the case without the skew. Thus, the skew only changes the off diagonal 2x2s. I can now extract R from M.

Ran into a hitch: Why aren't my A matrices in proper form? I have to rotate them! (Of course. I should have known this already.) I spent the night developing a 2x2 rotation algorithm on MatLab. This should be easily extendible to the 4x4 case.

Also spent some time learning MatLab MIA commands.

QUESTIONS
1. I'm still confused about orbits and 4 degrees of freedom. Given an orbit matrix Z (4 columns), I can find a one turn map, M, by M = Z(i+1)*Z(i)^-1.

2003-07-09T09:34:00.000-07:00

Test Post. Please Disregard.