Repository: atcollab/at (language: python)

Issue:
"atsummary", "ringpara" and "atx" use "atlinopt6" for optics
This PR updates `atsummary`, `ringpara` and `atx` to use `atlinopt6` for optics. As a consequence:
- the input lattice may be either 4d or 6d ("radiation" ON or OFF),
- the off-momentum is properly taken into account, including for 6d lattices.

`atx` has deeper modifications:
- Since the input ring may be directly entered as a 6d lattice, optional syntaxes for configuring the "radiation ON" conversion are now useless. The are kept for compatibility but are deprecated,
- The optics method can be specified as a `'method'` keyword, with the following possibilities:
  - `@atlinopt6` is now the default value, applied to the 4d version of the lattice,
  - `@atlinopt4` may be used to apply the same method as the old version (it's equivalent to `atlinopt`),
  - `@atlinopt2` is available, but has very little interest.
  
  Some output fields are identical for the 3 methods (`alpha`, `beta`, `ClosedOrbit`, `Dispersion`…) but others are specific to each method,
- By default the linear optics is computed with the 4d version of the lattice, for speed and compatibility with the old version. However, with the `@atlinopt6` method, it can be computed with the 6d version by adding the `'6d'` flag.

This addresses some problems reported in #463. However there are still some incoherencies with the emittance. When off-momentum, the radiation integrals are computed with the corresponding modified values of betas and dispersion. However the result is still not in agreement with `atx`. A simpler example appears with the momentum compaction factor. It may be computed either directly by tracking (`mcf` function), or from the I<sub>1</sub> radiation integral. `atsummary` displays both results, which are in perfect agreement on-momentum, but differ off-momentum, suggesting that the I<sub>1</sub> integral is wrong. The only ingredients are dispersion and bending radius, so there is very little to do. Replacing the radius &rho; with &rho;*(1+&delta;) does not help…

Any idea, someone?


Discussion / hints:
Dear @lfarv 

It seems almost fine.
I do not understand why for esrf.m lattice

atx does not compute per default emittance, bunch length and so on.
I need to force atradon to get the expected result.

r = esrf;
atx(r) KO
atx(atradon(r)) OK
...
modemittance =
   NaN   NaN   NaN
modcoupling =
   NaN
projemittance =
   NaN   NaN   NaN
projcoupling =
   NaN
dampingtime =
   1.0e+07 *
  -0.002362334449781   0.408790442714825  -4.299921165235869
dampingJ =
   4.021027761813939  -0.023236874967468   0.002209113153529
fs =
     2.165810245594138e+03
espread =
     0
blength =
     0


Hello @lnadolski 

For some reason, the lattice generated by `esrf()` has active cavities but no radiation. That's why you get crazy damping times and no equilibrium emittance.

In the previous `atx` version, you had to provide a lattice **without longitudinal motion**, to get the optics with `atlinopt`, and it was internally switched to "radiation ON" with default options (RF, dipoles, quadrupoles) for `ohmienvelope`.

In this new version, the same thing happens if `'is_6d'` is false, but alternatively the user may provide a lattice **with longitudinal motion** (`'is_6d'` true). In this case, it is switched to "radiation OFF" for linear optics, and taken as it is for `ohmienvelope`. The advantage is that the user may tune the longitudinal motion as he likes. For instance here with the `esrf()` lattice: because of the cavity, `'is_6d'` is true but there is no damping.

I think this is how `atx` should behave:
- either you enter a lattice without longitudinal motion, and 6d is internally turned on with default options,
- or you enter a lattice with longitudinal motion, and you are responsible for the options.
However, you example shows one incoherency: the synchronous phase is computed with the energy loss per turn: `eloss=atgetU0(ron);`.
By default `atgetU0` uses the I<sub>2</sub> integral, which ignore the 6d state and takes all dipoles. So for the `esrf.m` lattice, `eloss` and `synchrophase` are not coherent with the other results.

`atgetU0` has another method using tracking which would be preferable here for coherence.

I'll look at that!
`atx`, `atsummary` and `ringpara` now get the energy loss with `atgetU0(ring,'method','tracking')` when the input has longitudinal motion ON. This gives losses according to the radiative elements declared in the lattice.
Dear @lfarv 

Thanks for you reply.
I understand your statement.
At the same time, this is a change of habits for the users.
Up to now, atsummary, atx, etc were providing all parameters as soon as an active cavity where in the lattice.
The default lattice is in the control room and in most our simulations is a lattice with an active cavity for 6D tracking, 6D closed orbit, etc. but without any damping.

I would prefer to automatically activate radiation if not present in atsummary just for the computation of the radiation integrals, the bunch length etc.

Best regards,

Laurent.
> Up to now, atsummary, atx, etc were providing all parameters as soon as an active cavity where in the lattice.

**Problem**: it's just the opposite: up to now `atsummary`, `atx`, etc need as input a lattice **without any longitudinal motion**. They need a cavity in the lattice to get the RF voltage, but it has to be inactive (`IdentityPass`). The reason is that the optics computation uses `atlinopt`, which forbids any change of momentum. See the help in `atlinopt`:
```
%
%ATLINOPT only works for CONSTANT energy. So no PassMethod should:
%   1. change the longitudinal momentum dP (cavities, magnets with radiation)
%   2. have any time dependence (localized impedance, fast kickers etc)
%
```
The way it works until now is:
- `atsummary` and `ringpara` return analytical results: radiation integrals, longitudinal motion... So they don't need to turn anything on, neither cavities nor radiation. They work only with non-coupled lattices. And as I said earlier, it seems that radiation integrals are not correct off-momentum, even with the adequate optical functions
- `atx` computes optics with radiation and cavities OFF, and then turns them on to call `ohmienvelope`. It works in the general case.

This PR does not change anything in this case: if there is no radiation or cavity, it works exactly the same as before and gives the same results.

What this PR adds is the ability to call these functions **with radiation or cavities ON**. And in that case, the attempt is to give results corresponding to the properties of the input lattice: with or without damping.

**Context**: the problem of using `atlinopt` with radiation or cavities is long-standing issue. @simoneliuzzo may say a lot about that. The fundamental cause is the fact that the presence of radiation or cavities is not checked. And if you violate the rule, you get results, usually not so wrong and even possibly right in some conditions (on-momentum, without radiation and if the cavity is the 1st element in the lattice, the results are correct). The consequences of this bad practice are really problematic:
- Adding a check preventing the use of 6d motion in `atlinopt` and its derivatives would break many things, possibly even within AT. So it is postponed, probably for ever,
- trying to give correct results either by internally turning radiation and cavities OFF or using `atlinopt6` which accepts longitudinal motion changes the results, which also creates problems (like here),
- leaving things unchanged  allows users to use AT in a way which gives "almost correct", but in fact wrong, results. It is not acceptable in my opinion.

I may add that the wrong naming of `atradon`, `atradoff` `"radiation ON"` etc. added to the confusion by hiding their action on cavities. This should be better with the aliases recently introduced.

**So now**: what to do with this PR:
- forget it and leave things as they are: as I said before, I don't like the possibility of allowing users to get wrong results, even usually "slightly" wrong.
- turn radiation and cavities OFF as the 1st step of `atx`, `atsummary` etc. It may slightly change the results, and it cancels the purpose of the PR: allowing computation with 6D on,
- keep it, forcing users using bad practices to change their habits: in that case, the simplest is to initially define the lattice with 6d turned off: no radiation and inactive cavity. And turn 6D on for the very few computations involving longitudinal motion: 6D tracking, 6D closed orbit, 6D transfer matrix, `ohmienvelope`.
- Any other idea, someone?

I understand your point, so we can take some time to find an acceptable solution…



Indeed, I have said a lot in the past, so I am not so keen to revive this discussion. 

In general I think the user should be free to : 
1) chose to run simulations with any combination of cavity and radiation states
2) set energy shifts in the lattice via dpp or cavity Frequency or both.

I personally do not know how to make this work.

May be accepting the (slightly) wrong results informing the user with a warning ('Warning: backward compatibility mode. The values of... and ... may be incorrect with RF .. and cavity .... . Please consider using <commands> instead.'), could be a solution. The motivation to accept wrong results would be to keep backward compatibility. 

Hello @simoneliuzzo !

> chose to run simulations with any combination of cavity and radiation states

That's exactly what this PR tries. But it introduces changes pointed out by @lnadolski : a lattice with a cavity and no radiation gives 0 losses, 0 (or pi) synchronous phase, no equilibrium emittance and infinite damping times: it's correct, but it was no such before…


> set energy shifts in the lattice via dpp or cavity Frequency or both.

It's too bad you are leaving now: I hope you would appreciate the next PR to be proposed right after this one: specify everywhere the off-momentum through one of these 3 keywords: 
- dp=* (off-momentum)
- dct=* (path lengthening)
- df=* (RF frequency shift)

A single one, to avoid contradicting requirements. 
> So now: what to do with this PR:
>
> - forget it and leave things as they are: as I said before, I don't like the possibility of allowing users to get wrong results, even usually "slightly" wrong.
> - turn radiation and cavities OFF as the 1st step of `atx`, `atsummary` etc. It may slightly change the results, and it cancels the purpose of the PR: allowing computation with 6D on,
> - keep it, forcing users using bad practices to change their habits: in that case, the simplest is to initially define the lattice with 6d turned off: no radiation and inactive cavity. And turn 6D on for the very few computations involving longitudinal motion: 6D tracking, 6D closed orbit, 6D transfer matrix, `ohmienvelope`.
> - Any other idea, someone?

@lnadolski : what do we do with the PR? Is it acceptable for you as it is, or do we postpone or modify?

Fix patch:
diff --git a/atmat/atphysics/LinearOptics/atlinopt4.m b/atmat/atphysics/LinearOptics/atlinopt4.m
index 417e532a8..08c3d095d 100644
--- a/atmat/atphysics/LinearOptics/atlinopt4.m
+++ b/atmat/atphysics/LinearOptics/atlinopt4.m
@@ -57,6 +57,10 @@
 %   corresponding to a frequency shift. The resulting deltap/p is returned
 %   in the 5th component of the ClosedOrbit field.
 %
+% [...] = ATLINOPT4(...,'df',DF)
+%   Analyses the off-momentum lattice by specifying the RF frequency shift.
+%   The resulting deltap/p is returned in the 5th component of the ClosedOrbit field.
+%
 % [...] = ATLINOPT4(...,'orbit',ORBITIN)
 %   Do not search for closed orbit. Instead ORBITIN,a 6x1 vector
 %   of initial conditions is used: [x0; px0; y0; py0; DP; 0].
@@ -80,7 +84,7 @@
 %	[3] Brian W. Montague Report LEP Note 165, CERN, 1979
 %	[4] D.Sagan, D.Rubin Phys.Rev.Spec.Top.-Accelerators and beams, vol.2 (1999)
 %
-%  See also atlinopt atlinopt2 atlinopt6 tunechrom
+%  See also atlinopt2 atlinopt6 tunechrom
 
 % [...] = ATLINOPT4(...,'coupled',flag)
 %   Private keyword for computation without coupling (used for atlinopt2)
@@ -94,9 +98,9 @@
 [get_chrom,varargs]=getflag(varargin,'get_chrom');
 [get_w,varargs]=getflag(varargs,'get_w');
 [twiss_in,varargs]=getoption(varargs,'twiss_in',[]);
+[dp,varargs]=getoption(varargs,'dp',0.0);
+[dpargs,varargs]=getoption(varargs,{'orbit','dct','df'});
 [orbitin,varargs]=getoption(varargs,'orbit',[]);
-[ct,varargs]=getoption(varargs,'ct',NaN);
-[dp,varargs]=getoption(varargs,'dp',0);
 [coupled,varargs]=getoption(varargs,'coupled',true);
 [mkey,varargs]=getoption(varargs,'mkey','M');
 [DPStep,varargs]=getoption(varargs,'DPStep');
@@ -116,7 +120,7 @@
 end
 
 if isempty(twiss_in)        % Circular machine
-    [orbit,orbitin]=findorbit(ring,refpts,'dp',dp,'ct',ct,'orbit',orbitin,'XYStep',XYStep);
+    [orbit,orbitin]=findorbit4(ring,dp,refpts,dpargs{:},'XYStep',XYStep);
     dp=orbitin(5);
     [orbitP,o1P]=findorbit4(ring,dp+0.5*DPStep,refpts,orbitin,'XYStep',XYStep);
     [orbitM,o1M]=findorbit4(ring,dp-0.5*DPStep,refpts,orbitin,'XYStep',XYStep);
@@ -198,7 +202,7 @@
     'gamma',gamma, 'C',CL, 'A',AL, 'B',BL,...
     'beta',num2cell([BX,BY],2),...
     'alpha',num2cell([AX,AY],2),...
-    'mu',num2cell([MX,MY],2));
+    'mu',num2cell([MX,MY],2))';
 
 if get_w
     % Calculate optics for DP +/- DDP
diff --git a/atmat/atphysics/LinearOptics/atlinopt6.m b/atmat/atphysics/LinearOptics/atlinopt6.m
index 20803c1be..407330098 100644
--- a/atmat/atphysics/LinearOptics/atlinopt6.m
+++ b/atmat/atphysics/LinearOptics/atlinopt6.m
@@ -67,7 +67,10 @@
 %   Specify off-momentum 
 %
 % [...] = ATLINOPT6(...,'dct',DCT)
-%   Specify path lengthening
+%   Specify the path lengthening
+%
+% [...] = ATLINOPT6(...,'df',DF)
+%   Specify the RF frequency deviation
 %
 %  REFERENCES
 %   [1] Etienne Forest, Phys. Rev. E 58, 2481 – Published 1 August 1998
@@ -75,7 +78,7 @@
 %       Published 3 February 2006
 %   [3] Brian W. Montague Report LEP Note 165, CERN, 1979
 %
-%  See also atlinopt atlinopt2 atlinopt4 tunechrom
+%  See also atlinopt2 atlinopt4 tunechrom
 
 [ringdata,elemdata]=wrapper6d(ring,@xlinopt6,varargin{:});
 
@@ -83,7 +86,7 @@
         clight = PhysConstant.speed_of_light_in_vacuum.value;   % m/s
         [get_chrom, varargs]=getflag(varargin, 'get_chrom');
         [get_w, varargs]=getflag(varargs, 'get_w');
-        [dpargs,varargs]=getoption(varargs,{'orbit','dp','dct'});
+        [dpargs,varargs]=getoption(varargs,{'orbit','dp','dct','df'});
         [twiss_in,varargs]=getoption(varargs,'twiss_in',[]);
         [DPStep,~]=getoption(varargs,'DPStep');
         [cavargs,varargs]=getoption(varargs,{'cavpts'});
diff --git a/atmat/atphysics/LinearOptics/findm44.m b/atmat/atphysics/LinearOptics/findm44.m
index a55a6ddd9..56175ba10 100644
--- a/atmat/atphysics/LinearOptics/findm44.m
+++ b/atmat/atphysics/LinearOptics/findm44.m
@@ -1,4 +1,4 @@
-function [M44, varargout]  = findm44(LATTICE,dp,varargin)
+function [M44, varargout]  = findm44(LATTICE,varargin)
 %FINDM44 numerically finds the 4x4 transfer matrix of an accelerator lattice
 % for a particle with relative momentum deviation DP
 %
@@ -8,12 +8,12 @@
 %     (cavities , magnets with radiation, ...)
 %   2.have any time dependence (localized impedance, fast kickers, ...)
 %
-% M44 = FINDM44(LATTICE,DP) finds the full one-turn
+% M44 = FINDM44(LATTICE) finds the full one-turn
 %    matrix at the entrance of the first element
 %    !!! With this syntax FINDM44 assumes that the LATTICE
 %    is a ring and first finds the closed orbit
 %
-% [M44,T] = FINDM44(LATTICE,DP,REFPTS) also returns
+% [M44,T] = FINDM44(LATTICE,REFPTS) also returns
 %    4x4 transfer matrixes  between entrance of
 %    the first element and each element indexed by REFPTS.
 %    T is 4x4xlength(REFPTS) 3 dimensional array
@@ -32,7 +32,17 @@
 %    it is also the entrance of the first element
 %    after 1 turn: T(:,:,end) = M
 %
-% [...] = FINDM44(...,'orbit',ORBITIN)
+% [...] = FINDM44(RING,...,'dp',DP)
+% [...] = FINDM44(RING,DP,REFPTS,...)  (Deprecated syntax)
+%   Computes for the off-momentum DP
+%
+% [...] = FINDM44(RING,...,'dct',DCT)
+%   Computes for the path lenghening specified by CT.
+%
+% [...] = FINDM44(RING,...,'df',DF)
+%   Computes for a deviation of RF frequency DF
+%
+% [...] = FINDM44(RING,...,'orbit',ORBITIN)
 % [...] = FINDM44(RING,DP,REFPTS,ORBITIN)  (Deprecated syntax)
 %   Do not search for closed orbit. Instead ORBITIN,a 6x1 vector
 %   of initial conditions is used: [x0; px0; y0; py0; DP; 0].
@@ -44,11 +54,6 @@
 %   Same as above except that matrices returned in T are full 1-turn
 %   matrices at the entrance of each element indexed by REFPTS.
 %
-% [...] = FINDM44(...,'ct',CT)
-%   Instead of computing the linear optics for  the specified DP/P,
-%   computes for the path lenghening specified by CT.
-%   The DP argument is ignored.
-%
 % [M44,T,orbit] = FINDM44(...)
 %   In addition returns the closed orbit at the entrance of each element
 %   indexed by REFPTS.
@@ -62,10 +67,6 @@
 %       --------------------------------------
 %              2*delta
 %
-%    with optimal differentiation step delta given by !!!! DO LATER
-%    The relative error in the derivative computed this way
-%    is !!!!!!!!!!!!!!!!! DO LATER
-%    Reference: Numerical Recipes.
 
 if ~iscell(LATTICE)
     error('First argument must be a cell array');
@@ -74,7 +75,9 @@
 [fullflag,varargs]=getflag(varargin,'full');
 [XYStep,varargs]=getoption(varargs,'XYStep');
 [orbitin,varargs]=getoption(varargs,'orbit',[]);
-[ct,varargs]=getoption(varargs,'ct',NaN);
+varargs=getdparg(varargs);
+[dp,varargs]=getoption(varargs,'dp',0.0);
+[dpargs,varargs]=getoption(varargs,{'dct','df'});
 [refpts,orbitin,varargs]=getargs(varargs,[],orbitin,'check',@(x) ~(ischar(x) || isstring(x))); %#ok<ASGLU>
 
 if islogical(refpts)
@@ -90,10 +93,8 @@
         dp=orbitin(5);
     end
     orbitin = [orbitin(1:4);dp;0];
-elseif isnan(ct)
-    [~,orbitin] = findorbit4(LATTICE,dp);
 else
-    [~,orbitin]=findsyncorbit(LATTICE,ct);
+    [~,orbitin]=findorbit4(LATTICE,'dp',dp,dpargs{:});
 end
 
 refs=setelems(refpts,NE+1); % Add end-of-lattice
diff --git a/atmat/atphysics/LongitudinalDynamics/atsetcavity.m b/atmat/atphysics/LongitudinalDynamics/atsetcavity.m
index 169ca446e..05ab5d859 100644
--- a/atmat/atphysics/LongitudinalDynamics/atsetcavity.m
+++ b/atmat/atphysics/LongitudinalDynamics/atsetcavity.m
@@ -15,17 +15,20 @@
 %   Set the cavity frequency to the nominal value according to
 %   circumference and harmonic number
 %
-%NEWRING=ATSETCAVITY(RING,...,'Frequency','nominal','dp',dp)
+%NEWRING=ATSETCAVITY(RING,...,'Frequency','nominal','dp',DP)
 %   Set the cavity frequency to the nominal value for the specified dp
 %
-%NEWRING=ATSETCAVITY(RING,...,'Frequency','nominal','dct',dct)
+%NEWRING=ATSETCAVITY(RING,...,'Frequency','nominal','dct',DCT)
 %   Set the cavity frequency to the nominal value for the specified dct
 %
+%NEWRING=ATSETCAVITY(RING,...,'Frequency','nominal','df',DF)
+%   Set the cavity frequency to the nominal value + df
+%
 %NEWRING=ATSETCAVITY(RING,...,'Voltage',VOLTAGE,...)
 %   Set the total voltage (all cells) [V]
 %   The voltage of each cavity is VOLTAGE / N_CAVITIES / PERIODICITY
 %
-%NEWRING=ATSETCAVITY(RING,...,'HarmNumber',h,...)
+%NEWRING=ATSETCAVITY(RING,...,'HarmNumber',H,...)
 %   Set the harmonic number
 %
 %NEWRING=ATSETCAVITY(RING,...,'TimeLag',TIMELAG,...)
@@ -61,7 +64,7 @@
 %     on radflag says whether or not we want radiation on, which affects
 %     synchronous phase.
 %
-%  See also atSetCavityPhase, atsetRFcavity, atradon, atradoff, atgetU0
+%  See also atSetCavityPhase, atsetRFcavity, atenable_6d, atdisable_6d, atgetU0
 
 % Speed of light
 CLIGHT=PhysConstant.speed_of_light_in_vacuum.value;
@@ -74,6 +77,7 @@
 [timelag,varargs]=getoption(varargs, 'TimeLag', []);
 [dp,varargs]=getoption(varargs,'dp',NaN);
 [dct,varargs]=getoption(varargs,'dct',NaN);
+[df,varargs]=getoption(varargs,'df',NaN);
 
 if isempty(cavpts)
     [ncells,cell_h,beta0,cavpts]=atGetRingProperties(ring,'Periodicity','cell_harmnumber','beta','cavpts');
@@ -96,7 +100,10 @@
         lcell=findspos(ring,length(ring)+1);
         frev=beta0*CLIGHT/lcell;
         if (ischar(frequency) || isstring(frequency)) && strcmp(frequency, 'nominal')
-            if isfinite(dct)
+            hh=props_harmnumber(harmcell,cell_h);
+            if isfinite(df)
+                frev = frev + df/hh;
+            elseif isfinite(dct)
                 frev=frev * lcell/(lcell+dct);
             elseif isfinite(dp)
                 % Find the path lengthening for dp
@@ -106,7 +113,7 @@
                 dct=orbitout(6);
                 frev=frev * lcell/(lcell+dct);
             end
-            frequency = frev * props_harmnumber(harmcell,cell_h);
+            frequency = hh * frev;
         else
             harmcell=round(frequency/frev);
         end
@@ -138,9 +145,9 @@
     if radflag
         U0=atgetU0(ring);
         timelag= (lcell/(2*pi*harmcell))*asin(U0/vring);
-        ring=atradon(ring);  % set radiation on. nothing if radiation is already on
+        ring=atenable_6d(ring);  % set radiation on. nothing if radiation is already on
     else
-        ring=atradoff(ring,'RFCavityPass');  % set radiation off. nothing if radiation is already off
+        ring=atdisable_6d(ring,'cavipass','');  % set radiation off. nothing if radiation is already off
         timelag=0;
     end
     ring=atsetfieldvalues(ring, cavpts, 'TimeLag', timelag);
diff --git a/atmat/atphysics/Orbit/findorbit.m b/atmat/atphysics/Orbit/findorbit.m
index 62bad892e..70d06949b 100644
--- a/atmat/atphysics/Orbit/findorbit.m
+++ b/atmat/atphysics/Orbit/findorbit.m
@@ -1,59 +1,68 @@
-function [orbs,orbitin]  = findorbit(ring,varargin)
-%FINDORBIT find the closed orbit
-%
-%Depending on the lattice, FINDORBIT will:
-% - use findorbit6 if radiation is ON,
-% - use findsyncorbit if radiation is OFF and ct is specified,
-% - use findorbit4 otherwise
-%
-%[ORBIT,FIXEDPOINT]=FINDORBIT(RING,REFPTS)
-%
-% ORBIT:        6xNrefs, closed orbit at the selected reference points
-% FIXEDPOINT:   6x1 closed orbit at the entrance of the lattice
-%
-%[...]=FINDORBIT(RING,...,'dp',DP) Specify the momentum deviation when
-%   radiation is OFF (default: 0)
-%
-%[...]=FINDORBIT(RING,...,'dct',DCT) Specify the path lengthening when
-%   radiation is OFF (default: 0)
-%
-%[...]=FINDORBIT(RING,...,'orbit',ORBIT) Specify the orbit at the entrance
-%   of the ring, if known. FINDORBIT will then transfer it to the reference points
-%
-%[...]=FINDORBIT(RING,...,'guess',GUESS) The search will start with GUESS as
-%   initial conditions
-%
-%For other keywords, see the underlying functions.
-%
-% See also FINDORBIT4, FINDSYNCORBIT, FINDORBIT6.
-
-[orbitin,varargs]=getoption(varargin,'orbit',[]);
-[refpts,varargs]=getargs(varargs,[],'check',@(arg) isnumeric(arg) || islogical(arg));
-[dp,varargs]=getoption(varargs,'dp',NaN);
-[dct,varargs]=getoption(varargs,'dct',NaN);
-if isempty(orbitin)
-    if check_radiation(ring)    % Radiation ON: 6D orbit
-        if isfinite(dp) || isfinite(dct)
-            warning('AT:linopt','In 6D, "dp" and "dct" are ignored');
-        end
-        orbitin=xorbit_6(ring,varargs{:});
-    elseif isfinite(dct)
-        orbitin=xorbit_ct(ring,dct,varargs{:});
-    else
-        orbitin=xorbit_dp(ring,dp,varargs{:});
-    end
-    args={'KeepLattice'};
-else
-    args={};
-end
-
-if islogical(refpts)
-    refpts=find(refpts);
-end
-if isempty(refpts)
-    % return only the fixed point at the entrance of RING{1}
-    orbs=orbitin;
-else
-    orbs=linepass(ring,orbitin,refpts,args{:});
-end
-end
+function [orbs,orbitin]  = findorbit(ring,varargin)
+%FINDORBIT find the closed orbit
+%
+%Depending on the lattice, FINDORBIT will:
+% - use findorbit6 if radiation is ON,
+% - use findsyncorbit if radiation is OFF and ct is specified,
+% - use findorbit4 otherwise
+%
+%[ORBIT,FIXEDPOINT]=FINDORBIT(RING,REFPTS)
+%
+% ORBIT:        6xNrefs, closed orbit at the selected reference points
+% FIXEDPOINT:   6x1 closed orbit at the entrance of the lattice
+%
+%[...]=FINDORBIT(RING,...,'dp',DP) Specify the momentum deviation when
+%   radiation is OFF (default: 0)
+%
+%[...]=FINDORBIT(RING,...,'dct',DCT) Specify the path lengthening when
+%   radiation is OFF (default: 0)
+%
+%[...]=FINDORBIT(RING,...,'df',DF) Specify the RF frequency deviation
+%   radiation is OFF (default: 0)
+%
+%[...]=FINDORBIT(RING,...,'orbit',ORBIT) Specify the orbit at the entrance
+%   of the ring, if known. FINDORBIT will then transfer it to the reference points
+%
+%[...]=FINDORBIT(RING,...,'guess',GUESS) The search will start with GUESS as
+%   initial conditions
+%
+%For other keywords, see the underlying functions.
+%
+% See also FINDORBIT4, FINDSYNCORBIT, FINDORBIT6.
+
+[orbitin,varargs]=getoption(varargin,'orbit',[]);
+[refpts,varargs]=getargs(varargs,[],'check',@(arg) isnumeric(arg) || islogical(arg));
+[dp,varargs]=getoption(varargs,'dp',NaN);
+[dct,varargs]=getoption(varargs,'dct',NaN);
+[df,varargs]=getoption(varargs,'df',NaN);
+if isempty(orbitin)
+    if check_6d(ring)    % Radiation ON: 6D orbit
+        if isfinite(dp) || isfinite(dct) || isfinite(df)
+            warning('AT:linopt','In 6D, "dp", "dct" and "df" are ignored');
+        end
+        orbitin=xorbit_6(ring,varargs{:});
+    elseif isfinite(df)
+        [cell_l,cell_frev,cell_h]=atGetRingProperties(ring,'cell_length',...
+            'cell_harmnumber','cell_revolution_frequency');
+        dct=-cell_l*df/(cell_frev*cell_h+df);
+        orbitin=xorbit_ct(ring,dct,varargs{:});
+    elseif isfinite(dct)
+        orbitin=xorbit_ct(ring,dct,varargs{:});
+    else
+        orbitin=xorbit_dp(ring,dp,varargs{:});
+    end
+    args={'KeepLattice'};
+else
+    args={};
+end
+
+if islogical(refpts)
+    refpts=find(refpts);
+end
+if isempty(refpts)
+    % return only the fixed point at the entrance of RING{1}
+    orbs=orbitin;
+else
+    orbs=linepass(ring,orbitin,refpts,args{:});
+end
+end
diff --git a/atmat/atphysics/Orbit/findorbit4.m b/atmat/atphysics/Orbit/findorbit4.m
index e8af2457e..ad2cbc67b 100644
--- a/atmat/atphysics/Orbit/findorbit4.m
+++ b/atmat/atphysics/Orbit/findorbit4.m
@@ -1,76 +1,90 @@
-function [orb4,orbitin] = findorbit4(ring,dp,varargin)
-%FINDORBIT4 finds closed orbit in the 4-d transverse phase
-% space by numerically solving  for a fixed point of the one turn
-% map M calculated with LINEPASS
-%
-%         (X, PX, Y, PY, dP2, CT2 ) = M (X, PX, Y, PY, dP1, CT1)
-%
-%    under the CONSTANT MOMENTUM constraint, dP2 = dP1 = dP and
-%    there is NO constraint on the 6-th coordinate CT
-%
-% IMPORTANT!!! FINDORBIT4 imposes a constraint on dP and relaxes
-%    the constraint on the revolution frequency. A physical storage
-%    ring does exactly the opposite: the momentum deviation of a
-%    particle on the closed orbit settles at the value
-%    such that the revolution is synchronous with the RF cavity
-%
-%                 HarmNumber*Frev = Frf
-%
-%    To impose this artificial constraint in FINDORBIT4
-%    PassMethod used for any elemen SHOULD NOT
-%    1. change the longitudinal momentum dP (cavities , magnets with radiation)
-%    2. have any time dependence (localized impedance, fast kickers etc)
-%
-% FINDORBIT4(RING,dP) is 4x1 vector - fixed point at the
-%    entrance of the 1-st element of the RING (x,px,y,py)
-%
-% FINDORBIT4(RING,dP,REFPTS) is 4xLength(REFPTS)
-%   array of column vectors - fixed points (x,px,y,py)
-%   at the entrance of each element indexed by the REFPTS array.
-%   REFPTS is an array of increasing indexes that  select elements
-%   from the range 1 to length(RING)+1.
-%   See further explanation of REFPTS in the 'help' for FINDSPOS
-%
-% FINDORBIT4(RING,dP,REFPTS,GUESS)
-% FINDORBIT4(...,'guess',GUESS)     The search for the fixed point
-%   starts from initial condition GUESS. Otherwise the search starts from
-%   [0; 0; 0; 0; 0; 0]. GUESS must be a 6x1 vector.
-%
-% FINDORBIT4(...,'orbit',ORBIT)     Specify the orbit at the entrance
-%   of the ring, if known. FINDORBIT4 will then transfer it to the
-%   reference points. ORBIT must be a 6x1 vector.
-%
-% [ORBIT, FIXEDPOINT] = FINDORBIT4( ... )
-%	The optional second return parameter is a 6x1 vector:
-%   closed orbit at the entrance of the RING.
-%
-% See also FINDSYNCORBIT, FINDORBIT6.
-
-if ~iscell(ring)
-    error('First argument must be a cell array');
-end
-[orbitin,varargs]=getoption(varargin,'orbit',[]);
-[dct,varargs]=getoption(varargs,'dct',NaN);
-[refpts,varargs]=getargs(varargs,[],'check',@(arg) isnumeric(arg) || islogical(arg));
-if isempty(orbitin)
-    if isfinite(dct)
-        orbitin=xorbit_ct(ring,dct,varargs);
-    else
-        orbitin=xorbit_dp(ring,dp,varargs);
-    end
-    args={'KeepLattice'};
-else
-    args={};
-end
-
-if islogical(refpts)
-    refpts=find(refpts);
-end
-if isempty(refpts)
-    % return only the fixed point at the entrance of RING{1}
-    orb4=orbitin(1:4,1);
-else
-    orb6 = linepass(ring,orbitin,refpts,args{:});
-    orb4 = orb6(1:4,:);
-end
+function [orb4,orbitin] = findorbit4(ring,varargin)
+%FINDORBIT4 finds closed orbit in the 4-d transverse phase
+% space by numerically solving  for a fixed point of the one turn
+% map M calculated with LINEPASS
+%
+%         (X, PX, Y, PY, dP2, CT2 ) = M (X, PX, Y, PY, dP1, CT1)
+%
+%    under the CONSTANT MOMENTUM constraint, dP2 = dP1 = dP and
+%    there is NO constraint on the 6-th coordinate CT
+%
+% IMPORTANT!!! FINDORBIT4 imposes a constraint on dP and relaxes
+%    the constraint on the revolution frequency. A physical storage
+%    ring does exactly the opposite: the momentum deviation of a
+%    particle on the closed orbit settles at the value
+%    such that the revolution is synchronous with the RF cavity
+%
+%                 HarmNumber*Frev = Frf
+%
+%    To impose this artificial constraint in FINDORBIT4
+%    PassMethod used for any elemen SHOULD NOT
+%    1. change the longitudinal momentum dP (cavities , magnets with radiation)
+%    2. have any time dependence (localized impedance, fast kickers etc)
+%
+%FINDORBIT4(RING) is 4x1 vector - fixed point at the
+%    entrance of the 1-st element of the RING (x,px,y,py)
+%
+%FINDORBIT4(RING,REFPTS) is 4xLength(REFPTS)
+%   array of column vectors - fixed points (x,px,y,py)
+%   at the entrance of each element indexed by the REFPTS array.
+%   REFPTS is an array of increasing indexes that  select elements
+%   from the range 1 to length(RING)+1.
+%   See further explanation of REFPTS in the 'help' for FINDSPOS
+%
+%FINDORBIT4(RING,DP,REFPTS,...) Obsolete syntax
+%FINDORBIT4(RING,...,'dp',DP)   Specifies the off-momentum
+%
+%FINDORBIT4(RING,...,'dct',DCT) Specifies the path lengthening
+%
+%FINDORBIT4(RING,...,'df',DF) Specifies RF frequency deviation
+%
+%FINDORBIT4(RING,dP,REFPTS,GUESS)
+%FINDORBIT4(...,'guess',GUESS)     The search for the fixed point
+%   starts from initial condition GUESS. Otherwise the search starts from
+%   [0; 0; 0; 0; 0; 0]. GUESS must be a 6x1 vector.
+%
+%FINDORBIT4(...,'orbit',ORBIT)     Specify the orbit at the entrance
+%   of the ring, if known. FINDORBIT4 will then transfer it to the
+%   reference points. ORBIT must be a 6x1 vector.
+%
+%[ORBIT, FIXEDPOINT] = FINDORBIT4( ... )
+%	The optional second return parameter is a 6x1 vector:
+%   closed orbit at the entrance of the RING.
+%
+% See also FINDSYNCORBIT, FINDORBIT6.
+
+if ~iscell(ring)
+    error('First argument must be a cell array');
+end
+[orbitin,varargs]=getoption(varargin,'orbit',[]);
+[dp,varargs]=getdparg(varargs,0.0);
+[dct,varargs]=getoption(varargs,'dct',NaN);
+[df,varargs]=getoption(varargs,'df',NaN);
+[refpts,varargs]=getargs(varargs,[],'check',@(arg) isnumeric(arg) || islogical(arg));
+if isempty(orbitin)
+    if isfinite(df)
+        [cell_l,cell_frev,cell_h]=atGetRingProperties(ring,'cell_length',...
+            'cell_harmnumber','cell_revolution_frequency');
+        dct=-cell_l*df/(cell_frev*cell_h+df);
+        orbitin=xorbit_ct(ring,dct,varargs{:});
+    elseif isfinite(dct)
+        orbitin=xorbit_ct(ring,dct,varargs{:});
+    else
+        orbitin=xorbit_dp(ring,dp,varargs{:});
+    end
+    args={'KeepLattice'};
+else
+    args={};
+end
+
+if islogical(refpts)
+    refpts=find(refpts);
+end
+if isempty(refpts)
+    % return only the fixed point at the entrance of RING{1}
+    orb4=orbitin(1:4,1);
+else
+    orb6 = linepass(ring,orbitin,refpts,args{:});
+    orb4 = orb6(1:4,:);
+end
 end
\ No newline at end of file
diff --git a/atmat/atphysics/Orbit/findsyncorbit.m b/atmat/atphysics/Orbit/findsyncorbit.m
index e5d686526..f4686f410 100644
--- a/atmat/atphysics/Orbit/findsyncorbit.m
+++ b/atmat/atphysics/Orbit/findsyncorbit.m
@@ -1,75 +1,89 @@
-function [orb5,orbitin] = findsyncorbit(ring,dct,varargin)
-%FINDSYNCORBIT finds closed orbit, synchronous with the RF cavity
-% and momentum deviation dP (first 5 components of the phase space vector)
-% by numerically solving  for a fixed point
-% of the one turn map M calculated with LINEPASS
-%
-%       (X, PX, Y, PY, dP2, CT2 ) = M (X, PX, Y, PY, dP1, CT1)
-%
-%    under constraints CT2 - CT1 =  dCT = C(1/Frev - 1/Frev0) and dP2 = dP1 , where
-%    Frev0 = Frf0/HarmNumber is the design revolution frequency
-%    Frev  = (Frf0 + dFrf)/HarmNumber is the imposed revolution frequency
-%
-% IMPORTANT!!!  FINDSYNCORBIT imposes a constraint (CT2 - CT1) and
-%    dP2 = dP1 but no constraint on the value of dP1, dP2
-%    The algorithm assumes time-independent fixed-momentum ring
-%    to reduce the dimensionality of the problem.
-%    To impose this artificial constraint in FINDSYNCORBIT
-%    PassMethod used for any element SHOULD NOT
-%    1. change the longitudinal momentum dP (cavities , magnets with radiation)
-%    2. have any time dependence (localized impedance, fast kickers etc).
-%
-%
-% FINDSYNCORBIT(RING,DCT) is 5x1 vector - fixed point at the
-%		entrance of the 1-st element of the RING (x,px,y,py,dP)
-%
-% FINDSYNCORBIT(RING,DCT,REFPTS) is 5xLength(REFPTS)
-%   array of column vectors - fixed points (x,px,y,py,dP)
-%   at the entrance of each element indexed by the REFPTS array.
-%   REFPTS is an array of increasing indexes that  select elements
-%   from the range 1 to length(RING)+1.
-%   See further explanation of REFPTS in the 'help' for FINDSPOS
-%
-% FINDSYNCORBIT(RING,DCT,REFPTS,GUESS)
-% FINDSYNCORBIT(...,'guess',GUESS)     The search for the fixed point
-%   starts from initial condition GUESS. Otherwise the search starts from
-%   [0; 0; 0; 0; 0; 0]. GUESS must be a 6x1 vector.
-%
-% FINDSYNCORBIT(...,'orbit',ORBIT)     Specify the orbit at the entrance
-%   of the ring, if known. FINDSYNCORBIT will then transfer it to the
-%   reference points. ORBIT must be a 6x1 vector.
-%
-% [ORBIT, FIXEDPOINT] = FINDSYNCORBIT( ... )
-%	The optional second return parameter is a 6x1 vector:
-%   closed orbit at the entrance of the RING.
-%
-% See also FINDORBIT4, FINDORBIT6.
-%
-if ~iscell(ring)
-    error('First argument must be a cell array');
-end
-[orbitin,varargs]=getoption(varargin,'orbit',[]);
-[dp,varargs]=getoption(varargs,'dp',NaN);
-[refpts,varargs]=getargs(varargs,[],'check',@(arg) isnumeric(arg) || islogical(arg));
-if isempty(orbitin)
-    if isfinite(dp)
-        orbitin=xorbit_dp(ring,dp,varargs);
-    else
-        orbitin=xorbit_ct(ring,dct,varargs);
-    end
-    args={'KeepLattice'};
-else
-    args={};
-end
-
-if islogical(refpts)
-    refpts=find(refpts);
-end
-if isempty(refpts)
-    % return only the fixed point at the entrance of RING{1}
-    orb5=orbitin(1:5,1);
-else
-    orb6 = linepass(ring,orbitin,varargin{1},args{:});
-    orb5 = orb6(1:5,:);
-end
-end
+function [orb5,orbitin] = findsyncorbit(ring,varargin)
+%FINDSYNCORBIT finds closed orbit, synchronous with the RF cavity
+% and momentum deviation dP (first 5 components of the phase space vector)
+% by numerically solving  for a fixed point
+% of the one turn map M calculated with LINEPASS
+%
+%       (X, PX, Y, PY, dP2, CT2 ) = M (X, PX, Y, PY, dP1, CT1)
+%
+%    under constraints CT2 - CT1 =  dCT = C(1/Frev - 1/Frev0) and dP2 = dP1 , where
+%    Frev0 = Frf0/HarmNumber is the design revolution frequency
+%    Frev  = (Frf0 + dFrf)/HarmNumber is the imposed revolution frequency
+%
+% IMPORTANT!!!  FINDSYNCORBIT imposes a constraint (CT2 - CT1) and
+%    dP2 = dP1 but no constraint on the value of dP1, dP2
+%    The algorithm assumes time-independent fixed-momentum ring
+%    to reduce the dimensionality of the problem.
+%    To impose this artificial constraint in FINDSYNCORBIT
+%    PassMethod used for any element SHOULD NOT
+%    1. change the longitudinal momentum dP (cavities , magnets with radiation)
+%    2. have any time dependence (localized impedance, fast kickers etc).
+%
+%
+% FINDSYNCORBIT(RING) is 5x1 vector - fixed point at the
+%		entrance of the 1-st element of the RING (x,px,y,py,dP)
+%
+% FINDSYNCORBIT(RING,REFPTS) is 5xLength(REFPTS)
+%   array of column vectors - fixed points (x,px,y,py,dP)
+%   at the entrance of each element indexed by the REFPTS array.
+%   REFPTS is an array of increasing indexes that  select elements
+%   from the range 1 to length(RING)+1.
+%   See further explanation of REFPTS in the 'help' for FINDSPOS
+%
+%FINDORBIT4(RING,DCT,REFPTS,...) Obsolete syntax
+%FINDORBIT4(RING,...,'dct',DCT)  Specifies the path lengthening
+%
+%FINDORBIT4(RING,...,'dp',DP)   Specifies the off-momentum
+%
+%FINDORBIT4(RING,...,'df',DF) Specifies RF frequency deviation
+%
+% FINDSYNCORBIT(RING,DCT,REFPTS,GUESS)
+% FINDSYNCORBIT(...,'guess',GUESS)     The search for the fixed point
+%   starts from initial condition GUESS. Otherwise the search starts from
+%   [0; 0; 0; 0; 0; 0]. GUESS must be a 6x1 vector.
+%
+% FINDSYNCORBIT(...,'orbit',ORBIT)     Specify the orbit at the entrance
+%   of the ring, if known. FINDSYNCORBIT will then transfer it to the
+%   reference points. ORBIT must be a 6x1 vector.
+%
+% [ORBIT, FIXEDPOINT] = FINDSYNCORBIT( ... )
+%	The optional second return parameter is a 6x1 vector:
+%   closed orbit at the entrance of the RING.
+%
+% See also FINDORBIT4, FINDORBIT6.
+%
+if ~iscell(ring)
+    error('First argument must be a cell array');
+end
+[orbitin,varargs]=getoption(varargin,'orbit',[]);
+[dct,varargs]=getdparg(varargs,0.0,'key','dct');
+[dp,varargs]=getoption(varargs,'dp',NaN);
+[df,varargs]=getoption(varargs,'df',NaN);
+[refpts,varargs]=getargs(varargs,[],'check',@(arg) isnumeric(arg) || islogical(arg));
+if isempty(orbitin)
+    if isfinite(df)
+        [cell_l,cell_frev,cell_h]=atGetRingProperties(ring,'cell_length',...
+            'cell_harmnumber','cell_revolution_frequency');
+        dct=-cell_l*df/(cell_frev*cell_h+df);
+        orbitin=xorbit_ct(ring,dct,varargs{:});
+    elseif isfinite(dp)
+        orbitin=xorbit_dp(ring,dp,varargs{:});
+    else
+        orbitin=xorbit_ct(ring,dct,varargs{:});
+    end
+    args={'KeepLattice'};
+else
+    args={};
+end
+
+if islogical(refpts)
+    refpts=find(refpts);
+end
+if isempty(refpts)
+    % return only the fixed point at the entrance of RING{1}
+    orb5=orbitin(1:5,1);
+else
+    orb6 = linepass(ring,orbitin,refpts,args{:});
+    orb5 = orb6(1:5,:);
+end
+end
diff --git a/atmat/atphysics/TuneAndChromaticity/tunechrom.m b/atmat/atphysics/TuneAndChromaticity/tunechrom.m
index 45c2d835a..e2bf5e516 100644
--- a/atmat/atphysics/TuneAndChromaticity/tunechrom.m
+++ b/atmat/atphysics/TuneAndChromaticity/tunechrom.m
@@ -1,97 +1,100 @@
-function varargout = tunechrom(ring,varargin)
-%TUNECHROM computes linear tunes and chromaticities
-%
-%TUNE=TUNECHROM(RING)	Quick calculation of the fractional part of the tune
-%	from numerically computed transfer matrix.
-%
-% [TUNE, CHROM] = TUNECHROM(RINGD,DP,'get_chrom') - optionally computes the
-%    chromaticities by numerical differentiation from the difference between
-%   tune values at momentums DP+0.5*DPStep and DP-0.5*DPStep
-%
-%[...]=TUNECHROM(...,'orbit',ORBITIN	Do not search for closed orbit.
-%   Instead ORBITIN,a 6x1 vector of initial conditions is used:
-%   This syntax is useful to avoid recomputing the closed orbit if is
-%   already known;
-%
-%[...]=TUNECHROM(RING,DP)       (obsolete)
-%[...]=TUNECHROM(...,'dp',DP)	Specify the momentum deviation when
-%   radiation is OFF (default: 0)
-%
-%[...]=TUNECHROM(...,'dct',DCT) Specify the path lengthening when
-%   radiation is OFF (default: 0)
-%
-% Note: TUNECHROM computes tunes and chromaticities from the one-turn
-%   transfer matrix. The transfer matrix is computed from tracking using
-%   numerical differentiation. The error of numerical differentiation
-%   is sensitive to the step size. (Reference: Numerical Recipes)
-%   The calculation of tunes involves one numerical differentiation.
-%   The calculation of chromaticity involves TWO!!! numerical differentiations.
-%   The error in calculated chromaticity from may be substantial (~ 1e-5).
-%   Use the XYStep and DPStep keyword arguments to control the step size
-%   in chromaticity calculations
-%
-% See also ATLINOPT6
-[cpl1,allargs]=getflag(varargin,'coupling');	% ignored, kept for compatibility
-[cpl2,allargs]=getoption(allargs,'coupled',true);	% ignored, kept for compatibility
-allargs=getdparg(allargs);
-if cpl1 || ~cpl2
-    warning('AT:ObsoleteParameter','The "coupled" flag is ignored: coupling is always assumed');
-end
-[varargout{1:nargout}]=wrapper6d(ring,@xtunechrom,allargs{:});
-
-    function [tune, chrom] = xtunechrom(ring,is6d,varargin)
-        [oldchrom,varargs]=getflag(varargin,'chrom');
-        [chrom,varargs]=getflag(varargs,'get_chrom');
-        [dpargs,varargs]=getoption(varargs,{'orbit','dp','dct'});
-        [cavargs,varargs]=getoption(varargs,{'cavpts'});
-        [DPStep,~]=getoption(varargs,'DPStep');
-        if is6d
-            tunefunc=@tune6;
-        else
-            tunefunc=@tune4;
-        end
-
-        [~,orbitin]=findorbit(ring,dpargs{:},varargs{:});
-        
-        tune=tunefunc(ring,'orbit',orbitin,varargs{:});
-
-        if chrom || oldchrom || nargout == 2
-            if is6d
-                frf=get_rf_frequency(ring);
-                DFStep=-DPStep*mcf(atradoff(ring))*frf;
-                rgup=atsetcavity(ring,'Frequency',frf+0.5*DFStep,cavargs{:});
-                rgdn=atsetcavity(ring,'Frequency',frf-0.5*DFStep,cavargs{:});
-                [~,o1P]=findorbit(rgup,'guess',orbitin,varargs{:});
-                [~,o1M]=findorbit(rgdn,'guess',orbitin,varargs{:});
-                deltap=o1P(5)-o1M(5);
-            else
-                dp=orbitin(5);
-                [~,o1P]=findorbit4(ring,dp+0.5*DPStep,'guess',orbitin,varargs{:});
-                [~,o1M]=findorbit4(ring,dp-0.5*DPStep,'guess',orbitin,varargs{:});
-                deltap=DPStep;
-            end
-            tuneP=tunefunc(ring,'orbit',o1P,varargs{:});
-            tuneM=tunefunc(ring,'orbit',o1M,varargs{:});
-            chrom = (tuneP - tuneM)/deltap;
-        end
-
-        function f=get_rf_frequency(ring)
-            % Get the initial RF frequency
-            cavities=ring(atgetcells(ring, 'Frequency'));
-            freqs=atgetfieldvalues(cavities,'Frequency');
-            f=freqs(1);
-        end
-
-    function tune=tune4(ring,varargin)
-        mm=findm44(ring,NaN,varargin{:});
-        [~,vals]=amat(mm);
-        tune=mod(angle(vals)/2/pi,1);
-    end
-
-    function tune=tune6(ring,varargin)
-        mm=findm66(ring,varargin{:});
-        [~,vals]=amat(mm);
-        tune=mod(angle(vals)/2/pi,1);
-    end
-    end
-end
+function varargout = tunechrom(ring,varargin)
+%TUNECHROM computes linear tunes and chromaticities
+%
+%TUNE=TUNECHROM(RING)	Quick calculation of the fractional part of the tune
+%	from numerically computed transfer matrix.
+%
+% [TUNE, CHROM] = TUNECHROM(RINGD,DP,'get_chrom') - optionally computes the
+%    chromaticities by numerical differentiation from the difference between
+%   tune values at momentums DP+0.5*DPStep and DP-0.5*DPStep
+%
+%[...]=TUNECHROM(...,'orbit',ORBITIN	Do not search for closed orbit.
+%   Instead ORBITIN,a 6x1 vector of initial conditions is used:
+%   This syntax is useful to avoid recomputing the closed orbit if is
+%   already known;
+%
+%[...]=TUNECHROM(RING,DP)       (obsolete)
+%[...]=TUNECHROM(RING,...,'dp',DP)	Specify the momentum deviation when
+%   radiation is OFF (default: 0)
+%
+%[...]=TUNECHROM(RING,...,'dct',DCT) Specify the path lengthening when
+%   radiation is OFF (default: 0)
+%
+%[...]=TUNECHROM(RING,...,'df',DF) Specify the RF frequency deviation when
+%   radiation is OFF (default: 0)
+%
+% Note: TUNECHROM computes tunes and chromaticities from the one-turn
+%   transfer matrix. The transfer matrix is computed from tracking using
+%   numerical differentiation. The error of numerical differentiation
+%   is sensitive to the step size. (Reference: Numerical Recipes)
+%   The calculation of tunes involves one numerical differentiation.
+%   The calculation of chromaticity involves TWO!!! numerical differentiations.
+%   The error in calculated chromaticity from may be substantial (~ 1e-5).
+%   Use the XYStep and DPStep keyword arguments to control the step size
+%   in chromaticity calculations
+%
+% See also ATLINOPT6
+[cpl1,allargs]=getflag(varargin,'coupling');	% ignored, kept for compatibility
+[cpl2,allargs]=getoption(allargs,'coupled',true);	% ignored, kept for compatibility
+allargs=getdparg(allargs);
+if cpl1 || ~cpl2
+    warning('AT:ObsoleteParameter','The "coupled" flag is ignored: coupling is always assumed');
+end
+[varargout{1:nargout}]=wrapper6d(ring,@xtunechrom,allargs{:});
+
+    function [tune, chrom] = xtunechrom(ring,is6d,varargin)
+        [oldchrom,varargs]=getflag(varargin,'chrom');
+        [chrom,varargs]=getflag(varargs,'get_chrom');
+        [dpargs,varargs]=getoption(varargs,{'orbit','dp','dct','df'});
+        [cavargs,varargs]=getoption(varargs,{'cavpts'});
+        [DPStep,~]=getoption(varargs,'DPStep');
+        if is6d
+            tunefunc=@tune6;
+        else
+            tunefunc=@tune4;
+        end
+
+        [~,orbitin]=findorbit(ring,dpargs{:},varargs{:});
+        
+        tune=tunefunc(ring,'orbit',orbitin,varargs{:});
+
+        if chrom || oldchrom || nargout == 2
+            if is6d
+                frf=get_rf_frequency(ring);
+                DFStep=-DPStep*mcf(atradoff(ring))*frf;
+                rgup=atsetcavity(ring,'Frequency',frf+0.5*DFStep,cavargs{:});
+                rgdn=atsetcavity(ring,'Frequency',frf-0.5*DFStep,cavargs{:});
+                [~,o1P]=findorbit(rgup,'guess',orbitin,varargs{:});
+                [~,o1M]=findorbit(rgdn,'guess',orbitin,varargs{:});
+                deltap=o1P(5)-o1M(5);
+            else
+                dp=orbitin(5);
+                [~,o1P]=findorbit4(ring,dp+0.5*DPStep,'guess',orbitin,varargs{:});
+                [~,o1M]=findorbit4(ring,dp-0.5*DPStep,'guess',orbitin,varargs{:});
+                deltap=DPStep;
+            end
+            tuneP=tunefunc(ring,'orbit',o1P,varargs{:});
+            tuneM=tunefunc(ring,'orbit',o1M,varargs{:});
+            chrom = (tuneP - tuneM)/deltap;
+        end
+
+        function f=get_rf_frequency(ring)
+            % Get the initial RF frequency
+            cavities=ring(atgetcells(ring, 'Frequency'));
+            freqs=atgetfieldvalues(cavities,'Frequency');
+            f=freqs(1);
+        end
+
+    function tune=tune4(ring,varargin)
+        mm=findm44(ring,NaN,varargin{:});
+        [~,vals]=amat(mm);
+        tune=mod(angle(vals)/2/pi,1);
+    end
+
+    function tune=tune6(ring,varargin)
+        mm=findm66(ring,varargin{:});
+        [~,vals]=amat(mm);
+        tune=mod(angle(vals)/2/pi,1);
+    end
+    end
+end
diff --git a/atmat/atutils/getdparg.m b/atmat/atutils/getdparg.m
index 657495923..5c5821ffe 100644
--- a/atmat/atutils/getdparg.m
+++ b/atmat/atutils/getdparg.m
@@ -2,7 +2,7 @@
 %GETDPARG Handle positional dp arguments
 %
 %[DP,VARARGS]=GETDPARG(VARARGS)
-%   If the 1st argument in VARARGS is empty or scalar numeric, it is
+%   If the 1st argument in VARARGS is a scalar numeric less than 1, it is
 %   considered as DP and removed from VARARGS. 
 %
 %VARARGS=GETDPARG(VARARGS)
@@ -10,20 +10,24 @@
 %   {'DP', DP} is added to VARARGS
 
 [check,varargs]=getoption(varargin,'check',@isdparg);
-[dp,varargs]=getargs(varargs,NaN); %#ok<ASGLU> 
-[dp,args]=getoption(args,'dp',dp);
-[dp,args]=getargs(args,dp,'check',check);
-if isempty(dp), dp=NaN; end
+[key,varargs]=getoption(varargs,'key','dp');
+[value,~]=getargs(varargs,NaN);     % Default value
+if ~isempty(args) && check(args{1}) % Positional dp argument
+    value=args{1};
+    args=args(2:end);
+end
+[value,args]=getoption(args,key,value);  % Keyword argument takes precedence
 if nargout == 1
-    if isfinite(dp)
-        args=[args {'dp', dp}];
+    if isfinite(value)
+        args=[args {key, value}];
     end
     varargout={args};
 else
-    varargout={dp,args};
+    varargout={value,args};
 end
 
-    function v=isdparg(arg)
-        v=isempty(arg) || (isnumeric(arg) && isscalar(arg));
+    function ok=isdparg(arg)
+        % ~(arg >= 1) : NaN is accepted as a valid dp value
+        ok=isnumeric(arg) && isscalar(arg) && ~(arg >= 1);
     end
 end
\ No newline at end of file


Test patch:
diff --git a/atmat/attests/pytests.m b/atmat/attests/pytests.m
index 535af699b..2866f6776 100644
--- a/atmat/attests/pytests.m
+++ b/atmat/attests/pytests.m
@@ -101,7 +101,7 @@ function m44(testCase,lat2,dp)
             a2=cell(lattice.p.find_m44(dp));
             [pm44,~]=deal(a2{:});
             pm44=double(pm44);
-            testCase.verifyEqual(mm44,pm44,AbsTol=2.e-9);
+            testCase.verifyEqual(mm44,pm44,AbsTol=3.e-9);
         end
 
         function m66(testCase,lat2)
@@ -129,7 +129,7 @@ function tunechrom4(testCase,lat,dp)
             ptune=double(lattice.p.get_tune(pyargs(dp=dp)));
             pchrom=double(lattice.p.get_chrom(pyargs(dp=dp)));
             testCase.verifyEqual(mtune,ptune,AbsTol=2.e-9);
-            testCase.verifyEqual(mchrom,pchrom,AbsTol=2.e-4);
+            testCase.verifyEqual(mchrom,pchrom,AbsTol=3.e-4);
         end
 
         function tunechrom6(testCase,lat2,dp)
diff --git a/pyat/test_matlab/test_cmp_linopt6.py b/pyat/test_matlab/test_cmp_linopt6.py
index 5e8d816d1..700cb8645 100644
--- a/pyat/test_matlab/test_cmp_linopt6.py
+++ b/pyat/test_matlab/test_cmp_linopt6.py
@@ -62,7 +62,7 @@ def test_4d_analysis(engine, lattices, dp):
                  atol=5.e-4, rtol=0)
 
     _compare_6d(py_data, ml_data, fields, atol=1.e-6, rtol=1e-8)
-    _compare_6d(py_data, ml_data, [('W', 'W')], atol=1.e-6, rtol=1e-4)
+    _compare_6d(py_data, ml_data, [('W', 'W')], atol=1.e-6, rtol=2.e-4)
 
 
 @pytest.mark.parametrize('lattices',