hibtc / cpymad

Cython binding to MAD-X
http://hibtc.github.io/cpymad/
Other
27 stars 18 forks source link

cpymad--twiss problem #17

Closed Landau1908 closed 9 years ago

Landau1908 commented 9 years ago

Hi, Thomas,

I have two problems. 1/ The TWISS parameters can be calculated from a sequence named any names, for example 'x' as showed in the below. Why?

2/ If I want to plot a piece of lattice not the whole, such as range=qd[10]/qd[16], how I should do?

# -*- coding: utf-8 -*-
"""
Created on Wed Jun 10 14:22:58 2015

@author: Administrator
"""
from cpymad.madx import Madx

# create a new interpreter instance:
# the optional 'command_log' parameter can be used to store MAD-X
# command history.
madx = Madx(command_log="log.madx")

# determine the version of MAD-X that is actually loaded:
print(madx.version) 

# you execute arbitrary textual MAD-X commands:
#madx.input('call, file="input_file.madx";')

# there is a more convenient syntax available which does the same:
madx.command.call(file="spamatch_global.max")

# And for some commands there exist direct shortcuts:
#madx.call('/path/to/some/input_file.madx')

#Calculate TWISS parameters
twiss = madx.twiss(sequence='x')

#Your own analysis below
import matplotlib.pylab as plt
plt.plot(twiss['s'], twiss['betx'])
plt.plot(twiss['s'], twiss['bety'])
plt.plot(twiss['s'], twiss['Dx'])
plt.plot(twiss['s'], twiss['x'], 'g')
plt.show()
//sequence definition
// define the total length
circum=6912.0;

// define number of cells and therefore cell length
ncell = 108;
lcell = circum/ncell;

// forces and other constants;
// element definitions;

// define bending magnets and quadrupoles as multipoles

mbsps: multipole, knl={2.0*pi/(8*ncell)};

kqf =  4.36588E-02;
kqd = -4.36952E-02;
qf: multipole, knl:={0.0, kqf};
qd: multipole, knl:={0.0, kqd};

// define sextupoles for chromaticity correction
ksf =  2.04516700E-02;
ksd = -3.87252263E-02;
lsf:   multipole, knl:={0.0,0.0, ksf};
lsd:   multipole, knl:={0.0,0.0, ksd};

// define orbit correctors and beam position monitors
bpm: monitor, l=0.00001;
ch:  hkicker, l=0.00001;
cv:  vkicker, l=0.00001;

// sequence declaration;
cassps: sequence, refer=centre, l=circum;
start_machine: marker, at = 0;

// This defines ONE cell, repeat NCELL times
// to get the full machine
// SPS has 8 bending magnets per cell
   n = 1;
   while (n <= ncell) {
   qf: qf,   at=(n-1)*lcell+1.5425;;                   
   lsf: lsf,       at=(n-1)*lcell+3.6425;                
   ch:  ch,        at=(n-1)*lcell+4.2425;                
   bpm: bpm,       at=(n-1)*lcell+4.3425;                
   mbsps: mbsps,   at=(n-1)*lcell+5.0425;               
   mbsps: mbsps,   at=(n-1)*lcell+11.4425;                
   mbsps: mbsps,   at=(n-1)*lcell+23.6425;              
   mbsps: mbsps,   at=(n-1)*lcell+30.0425;                
   qd: qd,   at=(n-1)*lcell+33.5425;                  
   lsd: lsd,       at=(n-1)*lcell+35.6425;                    
   cv:  cv,        at=(n-1)*lcell+36.2425;                    
   bpm: bpm,       at=(n-1)*lcell+36.3425;                    
   mbsps: mbsps,   at=(n-1)*lcell+37.0425;                 
   mbsps: mbsps,   at=(n-1)*lcell+43.4425;                   
   mbsps: mbsps,   at=(n-1)*lcell+55.6425;                
   mbsps: mbsps,   at=(n-1)*lcell+62.0425;                  

   n = n + 1;
}
end_machine: marker at=circum;            
endsequence;
//madx definition
TITLE, s='MAD-X test';
call file="sps.seq";
option,-echo,-thin_foc;
// option,debug,-echo;

Beam, particle = proton, sequence=cassps, energy = 450.0,
          NPART=1.05E11, sige=      4.5e-4 ;

use, period=cassps;   

select,flag=twiss,column=name,s,x,y,mux,betx,muy,bety,dx,dy;
twiss,save,centre,file=twiss.out;
plot, haxis=s, vaxis=x, betx, bety,colour=100,
      range=qd[10]/qd[16];

match, sequence=cassps;   
   vary,name=kqf, step=0.00001;
   vary,name=kqd, step=0.00001;
   global,sequence=cassps,Q1=26.58;
   global,sequence=cassps,Q2=26.62;
   Lmdif, calls=10, tolerance=1.0e-21;
endmatch;
match, sequence=cassps;   
   vary,name=ksf, step=0.00001;
   vary,name=ksd, step=0.00001;
   global,sequence=cassps,DQ1=0.0;  
   global,sequence=cassps,DQ2=0.0;   
   Lmdif, calls=10, tolerance=1.0e-21;
endmatch;

select,flag=twiss,column=name,s,x,y,mux,betx,muy,bety,dx,dy;
twiss,save,centre,file=twiss.out;
plot, haxis=s, vaxis=x, betx, bety,colour=100,
      range=qd[10]/qd[16];

stop;
coldfix commented 9 years ago

Seriously, put more effort into it when reporting an issue!

Question 1

Short version: It doesn't work. Always make sure to use the correct sequence name! Read the MAD-X output carefully, it tells you that it doesn't work.

If you had created a minimal example you would have noticed this on your own (as soon as you remove the TWISS command from your madx file)!

Explanation: There is a TWISS command in your .madx file. MAD-X then stores the calculated values in a memory table named "twiss". When executing madx.twiss(sequence='x') from your python code, cpymad executes a command TWISS, sequence=x, .... For some reason MAD-X thinks that this is not an error and ignores the attempt. It does not signal the error so there is no way to infer this from python, and so the the python code just continues. It next loads values from the "twiss" table assuming its filled with happy new values. Since you previously calculated TWISS in some other way, there are indeed values, so this doesn't crash either. But it is an error!

Question 2

  1. you can use the range parameter of madx.twiss(). But it may not work when multiple elements have the same name.
  2. You can slice the array manually after figuring out what the index is. This looks like the following:
twiss = madx.twiss(...)
elements = madx.sequences['cassps'].elements
start = elements.index('qd:10')
stop = elements.index('qd:16')
twiss = {twiss[k][start:stop] for k in ['s', 'betx', 'bety']}
plt.plot(twiss['s'], twiss['betx'])
...

And by the way: why do you intermix the native madx twiss+plotting and cpymad? The whole purpose of cpymad is to access the madx in data in python, so you can use more sophisticated tools for data analysis, such as matplotlib.

Landau1908 commented 9 years ago

Very sorry for that.

Native madx twiss+plotting is just a model from madx guide, so I directly copy from that.

Could you give me an example about matching?

Landau1908 commented 9 years ago

Hi, Thomas

Your code has an error.

# -*- coding: utf-8 -*-
"""
Created on Wed Jun 10 14:22:58 2015

@author: Administrator
"""
from cpymad.madx import Madx

madx = Madx(command_log="log.madx")

# determine the version of MAD-X that is actually loaded:
print(madx.version) 

# there is a more convenient syntax available which does the same:
madx.command.call(file="spamatch_global.max")
#m = Madx(ElementList = 'cassps')
#Calculate TWISS parameters
twiss = madx.twiss(sequence='cassps')

elements = madx.sequences['cassps'].elements
start = elements.index('qd:10')
stop = elements.index('qd:16')
twiss = {twiss[k][start:stop] for k in ['s', 'betx', 'bety']}
#Your own analysis below
import matplotlib.pylab as plt
plt.plot(twiss['s'], twiss['betx'])
plt.plot(twiss['s'], twiss['bety'])
plt.plot(twiss['s'], twiss['Dx'])
plt.plot(twiss['s'], twiss['x'], 'g')
plt.plot(twiss['s'], twiss['betx'])
plt.show()

error information:

runfile('E:/cython_madx/test/cpymad/1.py', wdir='E:/cython_madx/test/cpymad')
MAD-X 5.02.04 (2014.11.14)
Traceback (most recent call last):

  File "<ipython-input-50-f4c790c83cb7>", line 1, in <module>
    runfile('E:/cython_madx/test/cpymad/1.py', wdir='E:/cython_madx/test/cpymad')

  File "C:\Python27\lib\site-packages\spyderlib\widgets\externalshell\sitecustomize.py", line 601, in runfile
    execfile(filename, namespace)

  File "C:\Python27\lib\site-packages\spyderlib\widgets\externalshell\sitecustomize.py", line 66, in execfile
    exec(compile(scripttext, filename, 'exec'), glob, loc)

  File "E:/cython_madx/test/cpymad/1.py", line 23, in <module>
    twiss = {twiss[k][start:stop] for k in ['s', 'betx', 'bety']}

  File "E:/cython_madx/test/cpymad/1.py", line 23, in <setcomp>
    twiss = {twiss[k][start:stop] for k in ['s', 'betx', 'bety']}

TypeError: unhashable type: 'numpy.ndarray'
coldfix commented 9 years ago

Sorry, the offending line must read:

twiss = {k: twiss[k][start:stop] for k in ['s', 'betx', 'bety']}
Landau1908 commented 9 years ago

Thanks. That works well. Haha

coldfix commented 9 years ago

There is a match function. Its docstring also contains a small example code, which can be accessed as follows:

help(madx.match)

However, its functionality is currently somewhat restricted. In particular, it doesn't support the global commands. You can either (a) contribute a patch to fix that, OR (b) write your match code like this

madx.command.match(sequence='cassps')
madx.command.vary(name='kqf', step=0.00001)
...
# 'global' is python keyword and therefore needs to be specified as string:
madx.command('global', sequence='cassps', Q1=26.58)
...
madx.command.endmatch()
Landau1908 commented 9 years ago

Can madx.command call any code in madx?

coldfix commented 9 years ago

Yes. As stated in the Usage.