1 | function [k, tab]=reson(ordre,period,fen) |
---|
2 | % function reson(ordre, period, fen) |
---|
3 | % Recherche des droites de resonance a x + b y = c |
---|
4 | % passant dans la fenetre [x1 x2 y1 y2] |
---|
5 | % Ordre de la resonance |a| + |b| |
---|
6 | % Critere de recherche c = period * entie |
---|
7 | % eg: |
---|
8 | % [k, tab] = reson(4,1,[18 19 10 11] |
---|
9 | % |
---|
10 | % Laurent S. Nadolski, Synchrotron SOLEIL, 02/04 |
---|
11 | |
---|
12 | %% check argument number |
---|
13 | if nargin < 2 |
---|
14 | help reson |
---|
15 | return |
---|
16 | end |
---|
17 | |
---|
18 | %% window dimensions |
---|
19 | x1=fen(1); x2=fen(2); y1=fen(3); y2=fen(4); |
---|
20 | |
---|
21 | tab = []; |
---|
22 | k=0; |
---|
23 | %% diagonal of the window used to check wether resonsonce is within the |
---|
24 | %% window frame |
---|
25 | d0 = sqrt(power(fen(2)-fen(1),2) + power(fen(4)-fen(3),2))/2; |
---|
26 | |
---|
27 | %% Look for resonance |
---|
28 | for m1 = -ordre:ordre |
---|
29 | am1 = abs(m1); |
---|
30 | for m2 = -ordre+am1:ordre-am1 |
---|
31 | %%% Recherche des bornes de C |
---|
32 | if (m2 == 0) |
---|
33 | c(1)=m1*x1; |
---|
34 | c(2)=m1*x2; |
---|
35 | elseif (m1/m2<=0) |
---|
36 | c(1)= m1*x1+m2*y2; |
---|
37 | c(2)= m1*x2+m2*y1; |
---|
38 | else |
---|
39 | c(1)= m1*x1+m2*y1; |
---|
40 | c(2)= m1*x2+m2*y2; |
---|
41 | end |
---|
42 | cmin = floor(min(c)); |
---|
43 | cmax= floor(max(c)); |
---|
44 | |
---|
45 | for m3=cmin:cmax |
---|
46 | if (mod(m3,period) == 0) |
---|
47 | if (abs(m1) + abs(m2) == ordre) |
---|
48 | % check whether resonance within the window |
---|
49 | d = distance(m1,m2,m3,fen); |
---|
50 | if (d < d0) |
---|
51 | k=k+1; |
---|
52 | tab(k,:)=[m1 m2 m3]; |
---|
53 | plot_reson(m1,m2,m3,fen) |
---|
54 | end |
---|
55 | |
---|
56 | end |
---|
57 | end |
---|
58 | end |
---|
59 | end |
---|
60 | end |
---|
61 | |
---|
62 | function d = distance(m1,m2,m3,fen) |
---|
63 | % calcul la distance entre la droite et la centre de la fenetre |
---|
64 | % centre de la fenetre |
---|
65 | pt= [fen(1) + (fen(2)-fen(1))/2, fen(3) + (fen(4)-fen(3))/2]; |
---|
66 | |
---|
67 | d = dist_point_droite(m1,m2,m3,pt); |
---|