source: MML/trunk/machine/SOLEIL/StorageRing/BBA/bba_last_sauvegarde_8sept07.m @ 4

Last change on this file since 4 was 4, checked in by zhangj, 11 years ago

Initial import--MML version from SOLEIL@2013
File size: 28.6 KB
Line 
1
2function [merit,CorConsEye,BpmRes,offs,OffsCor,center,ze,ord1,ord2,BpmEyeRes,CorDevListRes,iniCor,rmsFit,NameFile,NameFilePNG]=bba_last(RepRes,QuadFam,QuadDev,Plane,Method,nbDeltaCor,DeltaCor,offsetCor,retard,DeltaQp0,CorNumber)
3%% INPUTS
4%QuadFam = famille du quadrupole,
5%QuadDev = [cell element] qu quadrupole
6%Plane = 1 ou 2 ; horizontal ou vertical
7%Method = MEC (Most Effective Corrector) ou BMP4CM (Bump 4 cor),
8%nbDeltaCor = nb de pas du correcteur,
9%DeltaCor = Delta I d'un pas de correcteur,
10%offsetCor = offset en courant du correcteur pour recentrer la "parabole",
11%retard = délai en secondes pemettant d'attendre la fin de la montée en courant du correcteurs ou du quadrupole,
12%DeltaQp0 = Delta I pour un pas de quadrupole,
13%CorNumber = classement du correcteur en fonction de son efficacité (1=le plus efficace)
14
15%% EXEMPLE pour BBA M.E.C.
16% nbDeltaCor=13;
17% DeltaCor=0.1;% amps
18% offsetCor=0.7; %amps
19% retard=4;%20; sec
20% DeltaQp0=1.75;
21
22%% OUTPUTS
23% merit = tableau de 2 colonnes: colonne 1 -> courants dans le correcteur
24%                                colonne 2 -> valeur de la fonction de mérite
25%CorConsEye = relevé des valeur des consignes dans le correcteur -> dimension 1*nbDeltaCor
26%BpmRes = Valeur dans le BPM choisi (le plus proche)
27%OffsCor = offset à appliquer dans le correcteur pour centrer la parabole calculée
28%center = courant dans le correcteur pour lequel la parabole est minimale,
29%ze = abscisses calculées pour les fit (parabole et droite, 40 points),
30%ord1= ordonnées de la parabole évaluée en ze(i),
31%ord2= ordonnées de la droite évaluée en ze(i) -> pour le BPM choisi,
32%BpmEyeRes = [cell, element] du BPM choisi voir fonction "mec3",
33%CorDevListRes =  [cell, element] Correcteur choisi voir fonction "mec3",
34%iniCor = valeur initiale du correcteur
35%OffsetBpm = offset BPM en mm
36%rmsFit = variable permettant de qualifier le fit de la parabole: écart-type des distances entre le fit et les points mesurés
37%NameFile = nom du fichier dans lequel sont stockées toutes les données
38%NameFilePNG = nom du du fichier .PNG de la copie d'écran de l'interface graphique
39
40%BMP4COR
41DeltaPosMax=1;   %mm for BMP4CM Method
42
43
44condTune=0;      % permet la correction du point de fonctionnement QP pour revenir au nombre d'onde initial apres le BBA
45condQp=0;        % calcul automatique du DeltaQ -> 1 ou non -> 0 % pas disponible
46condBPMsigma=0;  % prise en compte ou non du sigma des Bpm -> 1 ou 0
47%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
48
49%switch2sim
50
51%setsp('HCOR',0)
52%setsp('VCOR',0)
53
54
55%Dx(s)=(-Dk(sQp) x(sQp))*[1/{1-k(sQp)*l Bet(sQp)/(2 tan(Pi nux))}]*
56%[rac(Bet(sBpm))*rac(Bet(sQp))]/(2 sin(pi nux))*
57%cos(fix(sBpm)-fix(sQp)-pi nux)
58% clear
59% famQ='Q1';
60% locQ=[4 1];
61% famB='BPMx';
62%
63% %position longitudinale du quadrupole considéré
64% sQp=getspos(famQ,locQ);
65% [famresB,locB,sBpm]=proche(famB,famQ,locQ);
66%
67% %extration du nombre d'onde
68% tune=gettune;
69% nux=tune(1);
70% nuz=tune(2);
71%
72% %extraction des fonctions bￜta au Qp considéré:
73% [bxQ,bzQ]=modelbeta(famQ,locQ);
74% %extraction des fonctions bￜta au BPM considéré:
75% [bxB,bzB]=modelbeta(famB,locB);
76%
77% %extraction des phases au Qp considéré:
78% [fixQ,fizQ]=modelphase(famQ,locQ);
79% %extraction des phases au BPM considéré:
80% [fixB,fizB]=modelphase(famB,locB);
81%% %%%%%%
82
83%switch2sim
84
85%setsp('HCOR',0)
86%setsp('VCOR',0)
87
88% % %% TEST DES VARIABLES D'ENTREES
89% % if ~iscellstr(QuadFam)
90% %     error('la famille entree n''est pas au format cell')
91% % end
92% %
93% % tailleQDev=size(QuadDev(1,:))
94% %
95% % if isempty(QuadDev)
96% %     error('Veuillez entrer le device du quadripole à étudier au format QuadDev=[cell elt]')
97% % elseif tailleQDev(2)~=2
98% %     error('Trop d''elements pour 1 device de QP:  QuadDev=[cell elt]')
99% % end
100% %
101% % if Plane<1 || Plane>2
102% %     error('Veuillez entrer le plan étudié: Horizontal: Plane=1 ; Vertical:Plane=2')
103% % end
104% %
105% % %ne fonctionne pas
106% % if ~ischar(Method)
107% %     error('Methodes disponibles au format string: ''MEC''')
108% % end
109% % %
110% %
111% %
112% % %% Mémorisation du répertoire de travail
113% % Rep    = getfamilydata('Directory','BBA');
114% %
115% % Annuncment = strcat('  ***         Registration Directory                             ***');
116% % disp('  ****************************************************************************');
117% % disp(Annuncment);
118% % disp('  **                                                                        **');
119% % disp('  ****************************************************************************');
120% % disp(' ')
121% % disp('  ****************************************************************************');
122% % disp('  **                             WARNING                                    **');
123% % disp('  **    New directory has to be defined after an offset application         **');
124% % disp('  **                                                                        **');
125% % disp('  ****************************************************************************');
126% %
127% % button = questdlg('Is the Registration Directory already existing ?','Registration Directory','yes','no','yes') ;
128% % if isequal(button,'yes')                            % la directory existe déjà
129% %     RepRes = uigetdir(Rep);                         % choix de la directory
130% % else
131% %     tmp2 = questdlg('Create directory ?','BBA','Yes','No','Yes');
132% %     if strcmpi(tmp2,'Yes')                          % create a new directory
133% %         prompt = {'Enter new directory name:'};
134% %         dlg_title = 'Input for Directory Name';
135% %         num_lines = 1;
136% %         def = {'work_december06'};
137% %         answer = inputdlg(prompt,dlg_title,num_lines,def);
138% %         RepRes = answer{1};
139% %         cd(Rep)                                    % on se place dans la directory amont /BBA
140% %         mkdir(RepRes)
141% %         RepRes = [Rep RepRes]                      % Nom complet de la directory
142% %     else                                           % alors on ne sait pas ce qu'on veut ??
143% %         disp('**  BBA cancelled  **')
144% %         return
145% %     end
146% % end
147% % disp('  ** BBA starting **')
148
149varfig=10; % utilité ??????
150tic;
151
152disp('%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%')
153%$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$%
154
155QuadFamilyList=QuadFam;
156QuadDevList=QuadDev;
157
158% Results are stored in a "dynamic' structure called BBA:
159% level1: BBA.(QuadFam)
160
161% Rappel du dernier fichier résultats de BBA
162[BBAnumb0, BBAnumbText0, NomText0]=SaveName(RepRes,0);
163
164
165%Chargement des fichiers nécessaire
166% A MODIFIER pour etre plus generique et prendre toujours la bonne matrice
167% reponse.
168load('respcor.mat') %réponse des correcteurs
169danow= datestr(now);
170texdat=danow(1:11);
171
172
173if BBAnumb0>0
174    %NameFile=[RepRes '/ResBBA_' BBAnumbText0 '.mat'];
175    NameFile=[RepRes '/' NomText0];  % dernier fichier BBA
176    %NameFile=[RepRes '/ResBBA_' BBAnumbText0 '_' texdat '.mat'];
177    load(NameFile);
178end
179%BBA;
180%définition du nom du nouveau fichier résultat BBA
181[BBAnumb1,BBAnumbText1,NomText1] = SaveName(RepRes,1);  %NomText1 n'est pas utulisé dans ce cas
182%NameFile=[RepRes '/ResBBA_' BBAnumbText1 '.mat'];
183NameFile=[RepRes '/ResBBA_' BBAnumbText1 '_' texdat '.mat'];  % nom du nouveau fichier BBA
184NameFilePNG=[RepRes '/SnapBBA_' BBAnumbText1 '_' texdat '.png'];
185%BBAnumb=BBAnumb+1;
186
187
188%[Reper '/DirBBA/resBBA_' BBAnumb '.mat']
189
190
191% recherche du tableau contenant l'écart type des BPM
192%rep='/home/matlabML/measdata/Ringdata/BPM'
193
194% NON UTILISE POUR L'INSTANT
195% rep = getfamilydata('Directory','BPMData')
196% col=dir(rep)
197%
198% for i=1:length(col)
199%     nom  = col(i).name;
200%     Lnom0= size(nom);
201%     Lnom = Lnom0(2);
202%     if  Lnom>30 & strcmpi(nom(Lnom-2:Lnom),'mat') & strcmpi(nom(1:7),'BPMData')
203%         nom2 = nom;
204%     end
205%
206% end
207% repBPM=char([rep '/' nom2]);
208% BPMdata=load(repBPM);
209%%%%%%%%%%%%%%%%%%%%%%%%%
210%%%%%%%%%%%%%%%%%%%%%%%%%
211
212
213
214%Plane=1;NameFilePNG
215%if      Plane==0
216%    loop=2;
217%    Plane=1;
218%else
219%    loop=1;
220%end
221%%
222%$$$$$$$$$$$
223%for m=1:loop
224%$$$$$$$$$$$
225
226
227
228if     Plane==1
229    BpmFam='BPMx';
230    CorFam='HCOR';
231    NamePlane='HPlane';
232    CorBpmResp=efficacy.HPlane; % matrice d'efficacité correcteur/BPM
233    %BpmMax=37; %non utilisé
234%    BpmSigma1=BPMdata.BPMxData.Sigma;
235elseif Plane==2
236    BpmFam='BPMz';
237    CorFam='VCOR';
238    NamePlane='VPlane';
239    CorBpmResp=efficacy.VPlane;
240    %BpmMax=4.5; %non utilisé
241   % BpmSigma1=BPMdata.BPMyData.Sigma;
242else
243    error('input must "1" for horizontal plane or "2" for vertical plane');
244end
245
246%if Method=='MEC'
247%% UTILISATION du MEC:
248%
249%mec(QuadFamilyList,QuadDev,Plane,CorBpmResp)
250%Demander si le fichier respcor.m convient
251%sinon [effic]=respcor(Plane)
252%$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$%
253% PREPARATION DES LISTES DE DEVICES %
254%$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$%
255
256BpmDevList = getlist(BpmFam);
257CorDevList = getlist(CorFam);
258
259for i=1:length(CorDevList)
260    CorDevListCell(i)=mat2cell(CorDevList(i,:),1,2);
261end
262
263
264
265%disp('poc');
266%getam('BPMx',[1 2])
267
268%$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$%
269% MEMORISATION DES VALEURS INITIALES %
270%$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$%
271
272% for i=1:length(BpmDevList) % voir peut ï¿œtre seulement Bpm Dispo???
273%     BpmDevList(i,:);
274%     InitValBpm(i)=getam(BpmFam,BpmDevList(i,:)); % MAT c'est long !! A modifier: pas trÚs performant
275%     BpmSigma0(i)=1;
276% end
277
278InitValBpm(1:length(BpmDevList))=getam(BpmFam,BpmDevList);
279BpmSigma0(1:120)=ones(1,120);
280
281if condBPMsigma==0
282    BpmSigma=BpmSigma0;
283else
284    BpmSigma=BpmSigma1;  % attention à définir
285end
286
287BpmSigma;
288InitValBpm;
289
290% for i=1:length(CorDevList) % voir peut ï¿œtre seulement COR Dispo???
291%     InitValCor(i)=getsp(CorFam,CorDevList(i,:)); % MAT c'est long ! A modifier: pas trÚs performant
292% end
293InitValCor(1:length(CorDevList))=getsp(CorFam,CorDevList); % MAT c'est long ! A modifier: pas trÚs performant
294
295
296% MODIFIER LA GESTION DES ERREURS: STATUS
297BPMStatus =family2status(BpmFam);
298CorStatus =family2status(CorFam);
299BpmFam;
300CorFam;
301
302if sum(BPMStatus)~=120
303    disp('Il y a des BPM qui ont un mauvais status')
304end
305
306if sum(CorStatus)~=56
307    disp('Il y a des COR qui ont un mauvais status')
308end
309%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
310
311
312
313
314taille=size(QuadFamilyList);
315
316disp('$$$$$$$$$$$$$$$$$$$$$$$$')
317for j=1:taille(1)
318
319    %% Choix de la variation minimale du courant dans le correcteur requise pour obtenir un bon BBA
320    QuadStatus=family2status(char(QuadFamilyList(j)),QuadDevList(j,:));
321
322    if QuadStatus~=1
323        disp('BBA impossible mauvais statut du QP')
324    end
325    InitValQp =getsp(char(QuadFamilyList(j)),QuadDevList(j,:));
326    %RECHERCHE DES ELEMENTS POUR LE BBA: BPM, CORRECTEURS
327    if strcmpi(Method,'MEC')
328        %CorNumb=1;
329        [CorFamRes,CorDevListRes,CorElemRes,EffRes,BpmFamRes,BpmEyeRes]=mec3(QuadFamilyList(j),QuadDevList(j,:),Plane,CorBpmResp,CorNumber);
330        CorStatus1 =family2status(CorFam,CorDevListRes);
331        BPMStatus1 =family2status(BpmFam,BpmEyeRes);
332
333        if CorStatus1~=1
334            disp('BBA impossible mauvais statut du COR')
335        end
336
337        if BPMStatus1~=1
338            disp('BBA impossible mauvais statut du BPM')
339        end
340    elseif strcmpi(Method,'BMP4CM')
341        [CorFamRes,CorDevListRes,BpmFamRes,BpmDevListRes,BpmEyeRes,ConsRes]=bmp4Cor(QuadFamilyList(j),QuadDevList(j,:),Plane);
342        tailleBp=size(BpmDevListRes);
343    end
344
345    tailleCor=size(CorDevListRes);
346    CorDevListRes;
347    for b=1:tailleCor(1)
348        InitValcorUsed(b)=getsp(CorFam,CorDevListRes(b,:)); % a modifier: peu performant
349    end
350    %%
351    %DETERMINATION DES TABLES DE VALEURS POUR LES COR,QP...
352    % deltaCor à ressortir de CorCurRange.m (pas encore programmé)
353    % nombre de deltaCor aussi à connaitre
354    clear CorCons CorSp BpmRes BpmVal QpCons QpSp
355    %nbDeltaCor=7;
356    %DeltaCor=0.75;%amps for MEC Method
357    % DeltaPosMax=1; %mm for BMP4CM Method
358    %retard=20;
359    var=0;
360    %stepsp(CorFam ,-nbDeltaCor*DeltaCor/2,CorDevListRes);
361    %stepsp(QuadFam,-nbDeltaQp*DeltaQp/4  ,QuadDev);
362    if strcmpi(Method,'MEC')
363        iniCor=InitValCor(CorElemRes);
364    elseif strcmpi(Method,'BMP4CM')
365        for i=1:4
366            inicorrect(i)=getam(CorFam,CorDevListRes(i,:));
367        end
368    end
369
370    for nbc=1:nbDeltaCor
371        if strcmpi(Method,'MEC')
372            CorTab(nbc)=InitValCor(CorElemRes)-((nbDeltaCor-1)*DeltaCor)/2+(nbc-1)*DeltaCor+offsetCor; %amps tableau ayant nbDeltaCor colonnes
373            % peut etre amélioré (sans boucle)
374        elseif strcmpi(Method,'BMP4CM')
375            %CorTab(nbc)=-DeltaPosMax/2+(nbc-1)*DeltaPosMax/(nbDeltaCor-1);
376
377            DeltaPos=DeltaPosMax/(nbDeltaCor-1);
378            CorTab(1)=-DeltaPosMax/2+offsetCor;
379            for valC=1:nbDeltaCor-1
380                CorTab(valC+1)=CorTab(1)+valC*DeltaPos;
381            end
382            %CorTab(3)=0.1;
383            % CorTab(4)=0.1;
384            %CorTab(5)=0.1;
385        end%mm   BMP4CM
386    end
387    CorTab;
388    % Calcul du deltaK pour obtenir une variation de DeltaPos
389    %DeltaQp=2;
390    nbDeltaQp=2;
391    % DeltaPos=1.5e-3; % A revoir
392    % DeltaQp=CalcDK(char(QuadFamilyList(j)),QuadDevList(j,:),BpmFamRes,BpmEyeRes,DeltaPos,Plane)*InitValQp
393
394    %% Choix du Delta I en fonction du DeltaNu MAX et du Delta I/I Max
395%     DeltaNu=0.01;BpmEyeRes
396%     DeltaI_I=0.01;
397%
398%     TuneTot=[-3.7e-3  9.37e-3 -2.33e-3   -8e-4  3.13e-3 -1.18e-3  6.75e-3 -1.88e-3 -7.95e-4  3.74e-3; ...
399%         2.91e-3 -2.74e-3  2.91e-3 3.74e-3   -1.68e-3  3.14e-3 -2.07e-3  2.14e-3  3.98e-3 -1.18e-3];
400%
401%     current=[125.72 161.32 75.63 147.71 202.58 118.45 213 180.55 179.35 210.48]*0.04;
402%     %[DK]=CalcDK(char(QuadFamilyList(j)),QuadDevList(j,:),BpmFamRes,BpmEyeRes,Plane)
403%
404%     for o=1:2
405%         for p=1:10
406%             TuneTot1(o,p)=DeltaNu/TuneTot(o,p);
407%             Dcourant(o,p)=min([abs(TuneTot1(o,p)) abs(current(p))]);
408%         end
409%     end
410%
411%
412%     numFamQ=char(QuadFamilyList(j));
413%     numQ=str2num(numFamQ(2:length(numFamQ)));
414%     DeltaQp=Dcourant(Plane,numQ);
415%     %DeltaQp=0.75;
416%     if condQp==1
417%         [DK]=CalcDK(char(QuadFamilyList(j)),QuadDevList(j,:),BpmFamRes,BpmEyeRes,Plane);
418%         DeltaQp=getsp(char(QuadFamilyList(j)),QuadDevList(j,:))*abs(DK);
419%     else
420        DeltaQp=DeltaQp0;
421%     end
422    %%
423
424    % Création des tableau de consignes pour les quadrupoles
425    % idée: toujours commencer "en haut" du cycle d'hystérésis
426    if InitValQp>0
427        QpTab(1)=InitValQp-DeltaQp;
428        QpTab(2)=InitValQp+DeltaQp;
429    elseif InitValQp<0
430        QpTab(1)=InitValQp+DeltaQp;
431        QpTab(2)=InitValQp-DeltaQp;
432    else
433        error('il y a un problÚme pour le calcul du DeltaK')
434    end
435 
436    %%
437    %
438
439    %sauvegarde des valeurs
440    NameQuadDev=char(['dev' num2str(QuadDevList(j,1)) '_'  num2str(QuadDevList(j,2))]);
441   
442    BBA.(char(QuadFamilyList(j))).(NameQuadDev).(char(NamePlane)).(Method).CorDev=CorDevListRes;
443    if strcmpi(Method,'MEC')
444        BBA.(char(QuadFamilyList(j))).(NameQuadDev).(char(NamePlane)).(Method).Eff   =EffRes;
445    end
446    BBA.(char(QuadFamilyList(j))).(NameQuadDev).(char(NamePlane)).(Method).BpmDev=BpmEyeRes;
447    BBA.(char(QuadFamilyList(j))).(NameQuadDev).(char(NamePlane)).(Method).Date  =datestr(now);
448
449
450    %sauvegarde des valeurs initiales
451    BBA.(char(QuadFamilyList(j))).(NameQuadDev).(char(NamePlane)).(Method).Init.Qp =InitValQp;
452    BBA.(char(QuadFamilyList(j))).(NameQuadDev).(char(NamePlane)).(Method).Init.Cor=InitValCor;
453    BBA.(char(QuadFamilyList(j))).(NameQuadDev).(char(NamePlane)).(Method).Init.Bpm=InitValBpm;
454
455    save(NameFile,'BBA');
456    load(NameFile);
457    %init=init+length(listfamQ);
458
459    toc;
460
461%% $$  Début du BBA   $$ %%
462    if strcmpi(Method,'MEC')
463
464        %%           DEBUT: boucle sur Qp (excursion +/- Dk)
465        vartext=0;
466        for v=1:nbDeltaQp
467            disp('Qp');
468            tuneNumber(v,:)=gettune;
469            current(v)=getdcct;
470            setsp(QuadFam,QpTab(v),QuadDev,-2);
471            pause(6);
472            %%   Debut: boucle sur COR (excursions nbDeltaCor)
473            clear CorCons CorSp
474            for u=1:nbDeltaCor
475                vartext=vartext+1;
476                %setsp(char(QuadFamilyList(j)),InitValQp,QuadDevList(j,:));
477                texto=char([num2str(vartext) '/' num2str(2*nbDeltaCor)]);
478                disp(texto);
479                drawnow;
480
481                setsp(CorFam,CorTab(u),CorDevListRes,-2);
482                pause(retard); % valeur définie par le programme en mode auto : valeur par défaut=2.5 sec
483                                                                %en mode manuel:valeur par défaut 4 sec
484                if u==1
485                    pause(9); % défini par l'expérience
486                end
487                %CorCons(u)=getsp(CorFam,CorDevListRes);
488                CorCons(u)=CorTab(u); % tentative 29 janvier
489                CorConsEye(u)=CorCons(u);
490                %CorSp(u)  =getsp(CorFam,CorDevListRes);
491
492                if v==1
493                    BpmRes(u)=getam(BpmFam,BpmEyeRes); %Mermorisation du BPM Eye
494                end
495
496                QpCons(u,v)=getsp(char(QuadFamilyList(j)),QuadDevList(j,:));
497                %QpSp(u,v)  =getsp(char(QuadFamilyList(j)),QuadDevList(j,:));
498                %tempo
499
500               
501                %Modifier, boucle inutile !!!
502                tic
503                BpmVal(1:120,v,u)=getam(BpmFam,BpmDevList);
504                %for w=1:length(BpmDevList)
505                    %BpmVal(w,v,u)=getam(BpmFam,BpmDevList(w,:)); %BpmVal(Bpm,Quad,Cor)
506                %end
507                toc
508                %%%%%%%%%%%%%%%%%%%%% test si l'orbite dépasse l'interlock en V  MAT 9 septembre 07
509%                 A = BpmVal(1:120,v,u) ;
510%                 I = find(abs(A)>0.7); % critÚre du test = 0.7 mm en absolu
511%                 if ~isempty(I)&Plane==2
512%                     disp(strcat('Quadrupole :',QuadFam,QuadDev)); % nom du quad et valeur du deltaQ
513%                     fprintf('delta I sur ce quad : %6.2f  A\n',DeltaQp);
514%                     disp( strcat('Correcteur :',QuadFam,QuadDev)); % nom du correcteur et valeur du deltacor
515%                     fprintf('delta I sur ce corr : %6.2f  A\n',DeltaCor);
516%                     fprintf('valeurs max de l''orbite %6.2f %6.2f mm\n',max(A(I)),min(A(I)));
517%                 end
518                %%%%%%%%%%%%%%%%%%%%%   
519               
520               
521
522            end
523            %%              FIN:boucle sur COR (excursions nbDeltaCor)
524            CorCons;
525            BpmRes;
526            QpCons;
527            %QpSp;
528            length(BpmVal);
529        end
530
531
532        %%
533
534    elseif strcmpi(Method,'BMP4CM')
535        iniCor=getam(BpmFam,BpmEyeRes);
536        %test du bump proposé
537        ocsi=setorbitbump(BpmFam,BpmDevListRes,ConsRes*CorTab(1),CorFam,[-2 -1 1 2],'Incremental','NoSetSP');
538        indCorRes=ocsi.CM.DeviceList; %doit etre le meme que CorDevListRes
539        test=sum(indCorRes-CorDevListRes);
540        if abs(test)~=0
541            error('bump impossible')
542        else
543            for u2=1:nbDeltaCor
544                ocsi=setorbitbump(BpmFam,BpmDevListRes,ConsRes*CorTab(u2),CorFam,[-2 -1 1 2],'Incremental','NoSetSP');
545                CorVal(:,u2)=ocsi.CM.Delta
546            end
547        end
548
549    vartext=0;
550        for v=1:nbDeltaQp
551            setsp(QuadFam,QpTab(v),QuadDev,-2);
552            tuneNumber(v,:)=gettune
553            current(v)=getdcct
554            pause(6);
555            disp('Qp');
556            %setsp(QuadFam,QpTab(v),QuadDev);
557            %%        Debut: boucle sur COR (excursions nbDeltaCor)
558            for u=1:nbDeltaCor
559                vartext=vartext+1;
560                texto=char([num2str(vartext) '/' num2str(2*nbDeltaCor)]);
561                disp(texto)
562                %setsp(char(QuadFamilyList(j)),InitValQp,QuadDevList(j,:));
563                for u1=1:4
564                    setsp(char(CorFam),CorVal(u1,u),CorDevListRes(u1,:),-2);
565                end
566                if u==1
567                    pause(6)
568                end
569                pause(retard);
570               
571                getsp('HCOR');
572
573                drawnow;
574                clear CorCons CorSp
575                %ocsi=setorbitbump(BpmFam,BpmDevListRes,ConsRes*CorTab(u),CorFam,[-2 -1 1 2],'Incremental','NoSetSP');
576                %indCorRes=ocsi.CM.DeviceList; %doit etre le meme que CorDevListRes
577                %test=sum(indCorRes-CorDevLisRes)
578                if abs(test)~=0
579                    error('bump impossible')
580                end
581                tailleCorRes=size(indCorRes);
582                for c=1:tailleCorRes(1)
583                    CorCons(u,c)=getsp(CorFam,indCorRes(c,:));
584                    %CorSp(u,c)  =getsp(CorFam,indCorRes(c,:));
585                end
586                if v==1
587                    BpmRes(u)=getam(BpmFam,BpmEyeRes); %Mermorisation du BPM Eye
588                end
589                QpCons(u,v)=getsp(char(QuadFamilyList(j)),QuadDevList(j,:));
590                %QpSp(u,v)  =getsp(char(QuadFamilyList(j)),QuadDevList(j,:));
591                %tempo
592
593                CorConsEye(u)=BpmRes(u);
594                tailleBList=size(BpmDevListRes);
595                ind1=dev2elem(BpmFam,BpmDevListRes(1,:));
596                ind2=dev2elem(BpmFam,tailleBList(1));
597                for w=1:ind1-1
598                    BpmVal(w,v,u)=getam(BpmFam,BpmDevList(w,:)); %BpmVal(Bpm,Quad,Cor)
599                end
600                for w=ind2+1:120
601                    BpmVal(w-tailleBList(1),v,u)=getam(BpmFam,BpmDevList(w,:));
602                end
603                %tempo
604            end
605        end
606    end
607
608
609
610    %%      Sauvegarde
611    BBA.(char(QuadFamilyList(j))).(NameQuadDev).(char(NamePlane)).(Method).Results.BpmData =BpmVal;
612    BBA.(char(QuadFamilyList(j))).(NameQuadDev).(char(NamePlane)).(Method).Results.BpmRes  =BpmRes;
613    BBA.(char(QuadFamilyList(j))).(NameQuadDev).(char(NamePlane)).(Method).Results.QuadCons=QpCons;
614    %BBA.(char(QuadFamilyList(j))).(NameQuadDev).(char(NamePlane)).(Method).Results.QuadSetp=QpSp;
615    BBA.(char(QuadFamilyList(j))).(NameQuadDev).(char(NamePlane)).(Method).Results.CorCons =CorCons;
616    BBA.(char(QuadFamilyList(j))).(NameQuadDev).(char(NamePlane)).(Method).Results.CorCons =CorCons;
617    BBA.(char(QuadFamilyList(j))).(NameQuadDev).(char(NamePlane)).(Method).Results.Current =current;
618    BBA.(char(QuadFamilyList(j))).(NameQuadDev).(char(NamePlane)).(Method).Results.Tune =tuneNumber;
619    %BBA.(char(QuadFamilyList(j))).(NameQuadDev).(char(NamePlane)).(Method).Results.CorSetp =CorSp;
620    save(NameFile,'BBA');
621    %%
622    %%      Retour aux Valeurs initiales Qp et COR
623    setsp(char(QuadFamilyList(j)),InitValQp,QuadDevList(j,:),-2); % Init Qp
624    for b=1:tailleCor(1)
625        setsp(CorFam,InitValcorUsed(b),CorDevListRes(b,:),-2);   % Init Cor
626    end
627    pause(retard);
628    tuneNumber(3,:)=gettune;
629    %%          Fin:boucle sur Qp (excursion +/- Dk)
630    if Plane==1
631        tuneWatch=tuneNumber(1,1);
632        tuneEnd  =tuneNumber(3,1);
633    elseif Plane==2
634        tuneWatch=tuneNumber(1,2);
635        tuneEnd  =tuneNumber(3,2);
636    end
637
638    %% boucle de correction du nombre d'onde
639    %Dnu0=tuneWatch-tuneEnd;
640    tuneEnd=gettune
641
642   
643    % par défaut non utilisé: condTune=0
644    % revoir pause
645    % revoir Delta I quad
646    if condTune==1
647        stepQp=signDQ*0.0005;
648        signDQ=sign(QpCons(1,1));
649        while abs(tuneWatch-tuneEnd)>1e-5
650            stepsp(char(QuadFamilyList(j)),stepQp,QuadDevList(j,:));
651            pause(0.5);
652            tuneEnd1=gettune;
653
654            if Plane==1
655                tuneEnd=tuneEnd1(1);
656            elseif Plane==2
657                tuneEnd=tuneEnd1(2);
658            end
659
660        end
661    end
662    disp('ole')
663
664    %%
665    %%      Recherche du minimum de la fonction de mérite
666    %load('resBBA.mat')
667    %data=BBA.(char(QuadFamilyList(j))).(NameQuadDev).(char(NamePlane)).MEC.Results.BpmData;
668    data=BpmVal;
669    datSize=size(data);
670    BpmSigma;
671
672    for r=1:datSize(3)     %itération sur les valeurs de correcteurs
673        point=0;
674        for t=1:datSize(1) %itération sur les valeurs des BPMs
675
676            point=point+1*((data(t,1,r)-data(t,2,r))/BpmSigma(t))^2;
677        end
678
679        merit(r,1)=CorConsEye(r);
680        merit(r,2)=1/datSize(1)*point;
681    end
682    BBA.(char(QuadFamilyList(j))).(NameQuadDev).(char(NamePlane)).(Method).Results.FitData =merit;
683    save(NameFile,'BBA');
684
685    %coef=polyfit(merit(:,1), merit(:,2),2); % calcul des coef de la parabole de la fonction de mï¿œrite
686    droite=polyfit(CorConsEye,BpmRes,1);        % calcul des coef de la droite du Bpm observé
687
688
689    %% Verif
690
691    %init:Moindres carrés
692 %%% A ALLEGER !!!!
693 
694    VectB=merit(:,2);
695    for q=1:length(merit)
696        MatFunc(q,1)=merit(q,1)^2;  MatFunc(q,2)=merit(q,1); MatFunc(q,3)=1;
697        M1(1,q)=MatFunc(q,1);% M1(1,2)=1; M1(1,3)=1;
698        M2(1,q)=MatFunc(q,2);%; M2(2,q)=1; M2(3,q)=0;
699        M3(1,q)=MatFunc(q,3);%; M3(2,q)=0; M3(3,q)=1;
700    end
701    M1;
702    ResVect(1,1)=M1*VectB;ResVect(2,1)=M2*VectB;ResVect(3,1)=M3*VectB;
703    ResMat(1,1)=M1*MatFunc*[1;0;0];ResMat(1,2)=M1*MatFunc*[0;1;0];ResMat(1,3)=M1*MatFunc*[0;0;1];
704    ResMat(2,1)=M2*MatFunc*[1;0;0];ResMat(2,2)=M2*MatFunc*[0;1;0];ResMat(2,3)=M2*MatFunc*[0;0;1];
705    ResMat(3,1)=M3*MatFunc*[1;0;0];ResMat(3,2)=M3*MatFunc*[0;1;0];ResMat(3,3)=M3*MatFunc*[0;0;1];
706    ResMat;
707    coef=inv(ResMat)*ResVect;
708
709    %%
710
711    center=-coef(2)/(2*coef(1));             % abscisse de l'offset BBA (amps)
712    offs=droite(1)*center+droite(2);         % ordonnï¿œe de l'offset BBA (mm)
713    %%      Fin de Recherche du minimum
714
715    BBA.(char(QuadFamilyList(j))).(NameQuadDev).(char(NamePlane)).(Method).Results.OffsetBpm =offs;
716    BBA.(char(QuadFamilyList(j))).(NameQuadDev).(char(NamePlane)).(Method).Results.CorCenter =center;
717
718    save(NameFile,'BBA');
719
720
721
722    for e=1:41
723        %ze(e)=min(BpmRes)+(e-1)*(max(BpmRes)-min(BpmRes))/40;
724        ze(e)=CorTab(1)+(e-1)*(CorTab(nbDeltaCor)-CorTab(1))/40;
725        %ze(e)=-2.5+(e-1)*5/40;
726        ord1(e)=parabole(ze(e),coef);
727        %ordco(e)=parabole(ze(e),valco);
728        ord2(e)=ligne(ze(e),droite);
729    end
730    toc;
731
732    for varTest=1:nbDeltaCor
733        ecart(varTest)=merit(varTest,2)-parabole(CorTab(varTest),coef);
734    end
735    rmsFit=qualFit(merit(:,1),merit(:,2),coef); % valeur affichée
736
737    BBA.(char(QuadFamilyList(j))).(NameQuadDev).(char(NamePlane)).(Method).Results.EcartFit =ecart;
738
739end
740
741
742
743
744%
745%Plane=2; % For running the prog on both Plane
746%%    %$$$$$$$$$$$$$$$
747%end
748NameFile;
749BpmRes;
750CorCons;
751%droite=polyfit(CorCons,BpmRes,1)
752iniCor;
753offs;
754
755ecart;
756CorDevListRes;
757if strcmpi(Method,'MEC')
758CorElemRes;
759EffRes;
760end
761disp('OffsCor');
762OffsCor=center-iniCor;
763
764
765% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
766% % Verif moindres carrés
767% xi=merit(:,1)
768% yi=merit(:,2)
769% xi0=0;
770% xi2=0;
771% xi3=0;
772% xi4=0;
773% mesur1=0;
774% mesur2=0;
775% mesur3=0;
776% for i=1:length(CorCons)
777%    xi4=xi4+xi(i)^4;
778%    xi3=xi3+xi(i)^3;
779%    xi2=xi2+xi(i)^2;
780%    xi0=xi0+1;
781%    mesur1=mesur1+yi(i)*xi(i)^2;
782%    mesur2=mesur2+yi(i)*xi(i);
783%    mesur3=mesur3+yi(i);
784% end
785%
786% mato(1,1)=  sum(xi4); mato(1,2)=sum(xi3);    mato(1,3)=  sum(xi2);
787% mato(2,1)=  sum(xi3); mato(2,2)=sum(xi2);    mato(2,3)=  sum(xi);
788% mato(3,1)=  sum(xi2); mato(3,2)=sum(xi);     mato(3,3)=  sum(xi0);
789%
790% mesur=[mesur1;mesur2;mesur3]/6
791%
792% mato=mato/6
793% valco=inv(mato)*mesur
794% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
795%
796%
797%
798%
799% for e=1:41
800%     ze(e)=CorTab(1)+(e-1)*(CorTab(nbDeltaCor)-CorTab(1))/40;
801%     %ze(e)=-2.5+(e-1)*5/40;
802%     ord1(e)=parabole(ze(e),coef);
803%     ordco(e)=parabole(ze(e),valco);
804%     ord2(e)=ligne(ze(e),droite);
805% end
806%
807% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
808%
809% y1=[min(BpmRes)-0.5 max(BpmRes)+0.5];
810% y2=[min(merit(:,2))-0.01 max(merit(:,2))+0.01];
811%
812% figure(1)
813% hold all
814% [Bx1,Hb1,Hb2]=plotyy(ze,ord2,ze,ord1,'plot')
815% %,'LineStyle','-.')
816% %textval=char([offs ' mm']);
817% %text(center,offs,textval,'HorizontalAlignment','left')
818% set(Hb1,'Color','black','LineStyle','-.');%,'LineStyle','--')
819% set(Hb2,'Color','black','LineStyle','-.');
820%
821%
822%
823%
824% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
825% %figure(2)
826%
827% [Ax1,Ha1,Ha2]=plotyy(CorCons,BpmRes,merit(:,1),merit(:,2) ,'plot');%,'Marker','x','Color','red')
828%
829% set(Ha1,'Marker','x','Color','blue');%,'LineStyle','--')
830% set(Ha2,'Marker','x','Color','red');
831%
832% x1=get(Ax1(1),'XLim')
833% set(Ax1,'Visible','off')
834% set(Ax1(1),'YLim',y1)
835% set(Ax1(2),'YLim',y2)
836%
837% pl1=plot([center center],[y1(1) offs])
838%
839% pl2=plot([x1(1) center],[offs offs])
840% pl3=plot(center,offs,'Marker','o','Color','black')%,text(center,offs,'\leftarrow ',...
841%      %'HorizontalAlignment','left'))
842%
843% set([pl1 pl2],'Color','black','LineStyle','-.');
844%
845%
846% set(Bx1,'Visible','on')
847% set(Bx1(1),'YLim',y1,'YColor','blue')
848% set(Bx1(2),'YLim',y2,'YColor','red')
849% xlabel('COR Value (amps)')
850% ylabel(Bx1(1),'(mm)')
851% ylabel(Bx1(2),'(merit)')
852%
853
854
855
856save('resu.mat','merit','ze','ord1','ord2','droite','coef','CorCons','BpmRes','coef','center','offs')
857
858%% Tout est enregistré dans une structure
859% BBA.(QuadFam).(['dev' QuadDev]).(Plane).(Method).CorDev
860%                                                 .eff
861%                                                 .BpmDev
862%                                                 .Date
863%                                                 .Results
864%...
865
866%example:
867% >> BBA.Q1.dev1_1.HPlane.MEC.CordeV
868% ans =
869% >> 5    1
Note: See TracBrowser for help on using the repository browser.