Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

So_fitting.f90 @ 430

Last change on this file since 430 was 430, checked in by touze, 11 years ago
import madx-5.01.00
File size: 95.7 KB

Line
1	!The Polymorphic Tracking Code
2	!Copyright (C) Etienne Forest and CERN
3
4	module S_fitting
5	USE MAD_LIKE
6	IMPLICIT NONE
7	public
8	PRIVATE FIND_ORBIT_LAYOUT, FIND_ORBIT_LAYOUT_noda
9	logical(lp), PRIVATE :: VERBOSE = .false.
10	integer :: max_fit_iter=20, ierror_fit=0, max_fiND_iter=40
11	real(dp) :: fuzzy_split=1.0_dp
12	real(dp) :: max_ds=0.0_dp
13	integer :: resplit_cutting = 0 ! 0 just magnets , 1 magnets as before / drifts separately
14	! 2 space charge algorithm
15	logical(lp) :: radiation_bend_split=my_false
16
17	INTERFACE FIND_ORBIT
18	! LINKED
19	! no use of TPSA
20	MODULE PROCEDURE FIND_ORBIT_LAYOUT_noda
21	END INTERFACE
22
23
24	contains
25	SUBROUTINE lattice_GET_CHROM(R,my_state,CHROM)
26	IMPLICIT NONE
27	TYPE(layout),target,INTENT(INOUT):: r
28	TYPE(internal_state), intent(in):: my_state
29	REAL(DP) CHROM(:)
30	TYPE(internal_state) state
31	real(dp) closed(6)
32	type(DAMAP) ID
33	TYPE(NORMALFORM) NORM
34	TYPE(REAL_8) Y(6)
35
36	STATE=((((my_state+nocavity0)+delta0)+only_4d0)-RADIATION0)
37
38	closed=0.0_dp
39	CALL FIND_ORBIT(R,CLOSED,1,STATE,1e-5_dp)
40	write(6,*) "closed orbit "
41	write(6,*) CLOSED
42	CALL INIT(STATE,2,0,BERZ)
43
44	CALL ALLOC(NORM)
45	CALL ALLOC(Y)
46	call alloc(id)
47
48	id=1
49	Y=CLOSED+id
50
51	CALL TRACK(R,Y,1,STATE)
52	NORM=Y
53	CHROM(1)=norm%DHDJ%V(1).SUB.'00001'
54	CHROM(2)=norm%DHDJ%V(2).SUB.'00001'
55	WRITE(6,*) "Fractional Tunes = ",norm%tune(1:2)
56	WRITE(6,*) "CHROMATICITIES = ",CHROM
57	CALL kill(NORM)
58	CALL kill(Y)
59	call kill(id)
60	if(.not.check_stable) then
61	CALL RESET_APERTURE_FLAG
62	write(6,*) " Flags were reset in lattice_GET_CHROM"
63	endif
64
65	end SUBROUTINE lattice_GET_CHROM
66
67
68	SUBROUTINE lattice_GET_tune(R,my_state,mf)
69	IMPLICIT NONE
70	TYPE(layout),target,INTENT(INOUT):: r
71	TYPE(internal_state), intent(in):: my_state
72	TYPE(internal_state) state
73	integer mf
74	real(dp) closed(6)
75	type(DAMAP) ID
76	TYPE(NORMALFORM) NORM
77	TYPE(REAL_8) Y(6)
78
79
80
81	STATE=my_state
82
83	closed=0.0_dp
84	CALL FIND_ORBIT(R,CLOSED,1,STATE,1e-5_dp)
85	write(6,*) "closed orbit "
86	WRITE(6,'(6(1x,g21.14))') CLOSED
87	CALL INIT(STATE,1,0,BERZ)
88
89	CALL ALLOC(NORM)
90	CALL ALLOC(Y)
91	call alloc(id)
92
93	id=1
94	Y=CLOSED+id
95
96	CALL TRACK(R,Y,1,STATE)
97	NORM=Y
98	closed=y
99	WRITE(6,'(6(1x,g21.14),a24)') CLOSED," <-- should be identical"
100	if(mf==6) then
101	WRITE(6,'(a19,3(1x,g21.14))') "Fractional Tunes = ",norm%tune(1:3)
102	if(norm%tune(3)/=0.0_dp.and.c_%ndpt==0) &
103	WRITE(6,'(a20,(1x,g21.14))') "Synchrotron period = ",1.d0/abs(norm%tune(3))
104	else
105	if(norm%tune(3)/=0.0_dp.and.c_%ndpt==0) then
106	WRITE(mf,'(4(1x,g21.14))') xsm0%ac%t/clight,norm%tune(1:3)
107	else
108	WRITE(mf,'(3(1x,g21.14))') xsm0%ac%t/clight,norm%tune(1:2)
109	endif
110	endif
111	CALL kill(NORM)
112	CALL kill(Y)
113	call kill(id)
114	if(.not.check_stable) then
115	CALL RESET_APERTURE_FLAG
116	write(6,*) " Flags were reset lattice_GET_tune"
117	endif
118
119	end SUBROUTINE lattice_GET_tune
120
121
122
123	SUBROUTINE compute_A_4d(r,my_state,filename,pos,del,no,MY_A)
124	IMPLICIT NONE
125	TYPE(layout),target,INTENT(INOUT):: r
126	TYPE(internal_state), intent(in):: my_state
127	integer pos,no,imod,ic,I
128	TYPE(internal_state) state
129	real(dp) closed(6),del
130	type(DAMAP) ID
131	TYPE(REAL_8) Y(6)
132	CHARACTER(*) FILENAME
133	TYPE(FIBRE), POINTER :: P
134	TYPE(DAMAP) MY_A
135	TYPE(NORMALFORM) NORM
136
137
138	STATE=((((my_state+nocavity0)-delta0)+only_4d0)-RADIATION0)
139
140	closed=0.0_dp
141	closed(5)=del
142	CALL FIND_ORBIT(R,CLOSED,pos,STATE,1e-5_dp)
143	write(6,*) "closed orbit "
144	write(6,*) CLOSED
145
146	CALL INIT(STATE,no,0,BERZ)
147
148	CALL ALLOC(Y); CALL ALLOC(MY_A); CALL ALLOC(NORM)
149	call alloc(id)
150	imod=r%n/10
151	id=1
152	Y=CLOSED+id
153	ic=0
154	p=>r%start
155	do i=1,r%n
156
157	CALL TRACK(R,Y,i,i+1,STATE)
158	if(mod(i,imod)==0) then
159	ic=ic+1
160	write(6,) ic10," % done "
161	endif
162	p=>p%next
163	enddo
164	id=y
165
166	call kanalnummer(ic)
167	open(unit=ic,file=FILENAME)
168
169
170	call print(ID,IC)
171	NORM=ID
172	MY_A=NORM%A_T
173	WRITE(6,*) " TUNES ",NORM%TUNE(1:2)
174	call print(MY_A,IC)
175	CLOSE(IC)
176
177
178	CALL kill(Y)
179	call kill(id)
180	call kill(NORM)
181	if(.not.check_stable) then
182	CALL RESET_APERTURE_FLAG
183	write(6,*) " Flags were reset in compute_A_4d"
184	endif
185
186	end SUBROUTINE compute_A_4d
187
188
189
190	SUBROUTINE compute_map_general(r,my_state,filename,pos,del,no)
191	IMPLICIT NONE
192	TYPE(layout),target,INTENT(INOUT):: r
193	TYPE(internal_state), intent(in):: my_state
194	TYPE(internal_state) state
195	integer pos,no,is,mf
196	real(dp) closed(6),del
197	type(DAMAP) ID
198	TYPE(REAL_8) Y(6)
199	CHARACTER(*) FILENAME
200	type(normalform) norm
201	type(taylor) betax,betax2
202	integer, allocatable :: expo1(:),expo2(:)
203	integer, allocatable :: expo3(:),expo4(:)
204
205	call kanalnummer(mf)
206	open(unit=mf,file=filename)
207
208	state=my_state !+nocavity0
209
210	if(state%nocavity) then
211	allocate(expo1(4),expo2(4))
212	allocate(expo3(4),expo4(4))
213	else
214	allocate(expo1(6),expo2(6))
215	allocate(expo3(6),expo4(6))
216	endif
217
218	expo1=0.0_dp;expo2=0.0_dp;
219	expo3=0.0_dp;expo4=0.0_dp;
220
221	closed=0.0_dp
222	closed(5)=del
223	CALL FIND_ORBIT(R,CLOSED,pos,my_state,1e-5_dp)
224	write(6,*) "closed orbit "
225	write(6,*) CLOSED
226
227	call init(state,no,c_%np_pol,berz)
228	call alloc(id); call alloc(norm);call alloc(y);call alloc(betax,betax2);
229
230	id=1; y=closed+id;
231	is=track_flag(r,y,pos,+state)
232
233	if(is/=0) then
234	write(6,*) "HELP"
235	stop
236	endif
237	id=y
238
239	call print(id,mf)
240	close(mf)
241	call kill(id); call kill(norm);call kill(y);call kill(betax,betax2);
242	deallocate(expo1,expo2)
243	deallocate(expo3,expo4)
244
245	END SUBROUTINE compute_map_general
246
247
248	SUBROUTINE compute_map_4d(r,my_state,filename,pos,del,no)
249	IMPLICIT NONE
250	TYPE(layout),target,INTENT(INOUT):: r
251	TYPE(internal_state), intent(in):: my_state
252	integer pos,no,imod,ic,I
253	TYPE(internal_state) state
254	real(dp) closed(6),del
255	type(DAMAP) ID
256	TYPE(REAL_8) Y(6)
257	CHARACTER(*) FILENAME
258	TYPE(FIBRE), POINTER :: P
259
260
261
262	STATE=((((my_state+nocavity0)-delta0)+only_4d0)-RADIATION0)
263
264	closed=0.0_dp
265	closed(5)=del
266	CALL FIND_ORBIT(R,CLOSED,pos,STATE,1e-5_dp)
267	write(6,*) "closed orbit "
268	write(6,*) CLOSED
269
270	CALL INIT(STATE,no,0,BERZ)
271
272	CALL ALLOC(Y)
273	call alloc(id)
274	imod=r%n/10
275	id=1
276	Y=CLOSED+id
277	ic=0
278	p=>r%start
279	do i=1,r%n
280
281	CALL TRACK(R,Y,i,i+1,STATE)
282	if(mod(i,imod)==0) then
283	ic=ic+1
284	write(6,) ic10," % done "
285	endif
286	p=>p%next
287	enddo
288	id=y
289
290	call kanalnummer(ic)
291	open(unit=ic,file=FILENAME)
292
293
294	call print(ID,IC)
295
296	CLOSE(IC)
297
298
299	CALL kill(Y)
300	call kill(id)
301
302	end SUBROUTINE compute_map_4d
303
304	SUBROUTINE FILL_BETA(r,my_state,pos,BETA,IB,DBETA,tune,tune2,a,ai,mat,clos)
305	IMPLICIT NONE
306	TYPE(layout),target,INTENT(INOUT):: r
307	TYPE(internal_state), intent(in):: my_state
308	REAL(DP), pointer :: BETA(:,:,:)
309	REAL(DP)DBETA,tune(:),tune2(:)
310	type(fibre),pointer :: p
311	integer i,IB,pos,mf
312	TYPE(internal_state) state
313	real(dp) closed(6)
314	type(DAMAP) ID
315	TYPE(NORMALFORM) NORM
316	TYPE(REAL_8) Y(6)
317	real(dp) dbetamax,db1,db2,db31,db32,eta1,eta2
318	real(dp),optional :: a(6,6),ai(6,6),mat(6,6),clos(6)
319
320	if(.not.associated(beta)) ALLOCATE(BETA(2,3,R%N))
321
322	eta1=0.0_dp
323	eta2=0.0_dp
324	STATE=((((my_state+nocavity0)-delta0)+only_4d0)-RADIATION0)
325	state=state+delta0
326
327	if(present(clos)) then
328	closed=clos
329	else
330	closed=0.0_dp
331	endif
332	if(pos/=0) then
333	CALL FIND_ORBIT(R,CLOSED,pos,STATE,1e-5_dp)
334	write(6,*) "closed orbit "
335	write(6,*) CLOSED
336	else
337	write(6,*) " Using a map "
338	endif
339	DBETA=0.0_dp
340	dbetamax=0.0_dp
341	CALL INIT(STATE,1,0,BERZ)
342
343	CALL ALLOC(NORM)
344	CALL ALLOC(Y)
345	call alloc(id)
346	if(pos/=0) then
347
348	id=1
349	Y=CLOSED+id
350
351	CALL TRACK(R,Y,1,STATE)
352	write(6,*) " stability ", c_%check_stable
353	id=y
354	else
355	call kanalnummer(mf)
356	open(unit=mf,file='map.dat')
357	call read(id,mf)
358	close(mf)
359	endif
360	NORM=id
361	if(present(a)) then
362	a=0.0_dp
363	ai=0.0_dp
364	a=norm%a_t
365	ai=norm%a_t**(-1)
366	id=y
367	mat=id
368	if(pos==0) Write(6,*) " Tunes ",norm%tune(1:2)
369	endif
370	if(pos/=0) then
371	if(ib==1) then
372	tune(1:2)=norm%tune(1:2)
373	endif
374	tune2(1:2)=norm%tune(1:2)
375	Write(6,*) " Tunes ", norm%tune(1:2)
376
377	y=closed+norm%a_t
378	p=>r%start
379	do i=pos,pos+r%n-1
380
381	CALL TRACK(R,Y,i,i+1,STATE)
382	beta(IB,1,i-pos+1)=(y(1).sub.'1')2 + (y(1).sub.'01')2
383	beta(iB,2,i-pos+1)=(y(3).sub.'001')2 + (y(3).sub.'0001')2
384	beta(IB,3,i-pos+1)=y(1).sub.'00001'
385
386	IF(IB==2) THEN
387	db1=ABS(beta(2,1,i-pos+1)-beta(1,1,i-pos+1))/beta(1,1,i-pos+1)
388	db2=ABS(beta(2,2,i-pos+1)-beta(1,2,i-pos+1))/beta(1,2,i-pos+1)
389	db31=ABS(beta(1,3,i-pos+1))
390	db32=ABS(beta(2,3,i-pos+1)-beta(1,3,i-pos+1))
391	DBETA=(db1+db2)/2.0_dp+dbeta
392	eta1=db31+eta1
393	eta2=db32+eta2
394	if( db1>dbetamax) dbetamax=db1
395	if( db2>dbetamax) dbetamax=db2
396	ENDIF
397	p=>p%next
398	enddo
399	DBETA=DBETA/R%N
400
401	IF(IB==2) WRITE(6,*) "<DBETA/BETA> = ",DBETA
402	IF(IB==2) WRITE(6,*) "MAXIMUM OF DBETA/BETA = ",dbetamax
403	IF(IB==2) WRITE(6,*) "<DETA/ETA> = ",eta2/eta1
404
405	endif
406	CALL kill(NORM)
407	CALL kill(Y)
408	call kill(id)
409	if(present(clos)) clos=closed
410
411	end SUBROUTINE FILL_BETA
412
413	SUBROUTINE comp_longitudinal_accel(r,my_state,no,h,filename)
414	IMPLICIT NONE
415	TYPE(layout),target,INTENT(INOUT):: r
416	TYPE(internal_state) state
417	real(dp) closed(6),l1,l2,p0c
418	type(DAMAP) ID
419	TYPE(NORMALFORM) NORM
420	TYPE(REAL_8) Y(6)
421	integer no,mf,i,h
422	CHARACTER(*) FILENAME
423	integer, allocatable :: j(:)
424	TYPE(internal_state), intent(in):: my_state
425	TYPE(fibre), pointer :: p
426	logical first
427	type(work) w
428	p=> r%start
429
430	l1=0.d0
431	l2=0.d0
432	first=.true.
433	w=p
434	p0c=w%p0c
435
436	call kanalnummer(mf,filename)
437
438	write(6,*)h,w%mass,w%p0c,w%beta0
439	write(mf,*) h
440	write(mf,*) w%mass,w%p0c,w%beta0
441
442	do i=1,r%n
443	if(p%mag%kind/=kind4) then
444	if(first) then
445	l1=l1+p%mag%p%ld
446	else
447	l2=l2+p%mag%p%ld
448	endif
449	else
450	if(.not.first) stop 111
451	l1=l1+p%mag%p%ld/2
452	l2=p%mag%p%ld/2
453	first=.false.
454	endif
455	p=>p%next
456	enddo
457	write(6,*) l1,l2
458	write(mf,*) l1,l2
459
460
461	STATE=my_state+nocavity0
462
463	CALL FIND_ORBIT(R,CLOSED,1,STATE,1e-5_dp)
464
465	CALL INIT(STATE,no,0)
466
467	CALL ALLOC(NORM)
468	CALL ALLOC(Y)
469	call alloc(id)
470
471	id=1
472	Y=CLOSED+id
473
474	CALL TRACK(R,Y,1,STATE)
475	NORM=Y
476
477	allocate(j(c_%nv))
478	j=0
479	do i=1,no
480	j(5)=i
481	write(6,*) norm%dhdj%v(3).sub.j
482	write(mf,*) norm%dhdj%v(3).sub.j
483	enddo
484
485	deallocate(j)
486
487	CALL kill(NORM)
488	CALL kill(Y)
489	call kill(id)
490	close(mf)
491
492	end SUBROUTINE comp_longitudinal_accel
493
494
495	SUBROUTINE comp_linear2(r,my_state,a,ai,mat,closed)
496	IMPLICIT NONE
497	TYPE(layout),target,INTENT(INOUT):: r
498	TYPE(internal_state), intent(in):: my_state
499	TYPE(internal_state) state
500	real(dp) closed(6),a(6,6),ai(6,6),mat(6,6)
501	type(DAMAP) ID
502	TYPE(NORMALFORM) NORM
503	TYPE(REAL_8) Y(6)
504
505
506	STATE=((((my_state+nocavity0)-delta0)+only_4d0)-RADIATION0)
507
508	CALL FIND_ORBIT(R,CLOSED,1,STATE,1e-5_dp)
509
510	CALL INIT(STATE,1,0,BERZ)
511
512	CALL ALLOC(NORM)
513	CALL ALLOC(Y)
514	call alloc(id)
515
516	id=1
517	Y=CLOSED+id
518
519	CALL TRACK(R,Y,1,STATE)
520	NORM=Y
521
522	Write(6,*) " Tunes ",norm%tune(1:2)
523	A=0.0_dp
524	AI=0.0_dp
525	MAT=0.0_dp
526	a=norm%a_t
527	AI=norm%a_t**(-1)
528	id=y
529	mat=id
530	CALL kill(NORM)
531	CALL kill(Y)
532	call kill(id)
533
534	end SUBROUTINE comp_linear2
535
536
537	subroutine lattice_fit_TUNE_gmap_auto(R,my_state,EPSF,TARG,name)
538	IMPLICIT NONE
539	TYPE(layout), target,intent(inout):: R
540	real(dp) , intent(IN),dimension(:)::TARG
541	real(dp) CLOSED(6)
542	TYPE(INTERNAL_STATE), intent(IN):: my_STATE
543	TYPE(INTERNAL_STATE) STATE
544	INTEGER I,SCRATCHFILE, MF
545	TYPE(TAYLOR), allocatable:: EQ(:)
546	TYPE(REAL_8) Y(6)
547	TYPE(NORMALFORM) NORM
548	integer :: neq=2, no=2,nt,j,it,NP
549	type(damap) id
550	type(gmap) g
551	TYPE(TAYLOR)t
552	character(nlp) name(2)
553	real(dp) epsf,epsr,epsnow,gam(2)
554	type(fibre), pointer:: p
555	logical dname(2)
556	TYPE(POL_BLOCK) poly(2)
557	call context(name(1))
558	call context(name(2))
559	poly(1)=0
560	poly(2)=0
561	poly(1)%ibn(2)=1
562	poly(2)%ibn(2)=2
563	! EPSF=.0001
564	SET_TPSAFIT=.FALSE.
565	dname=.false.
566	p=>r%start
567	do i=1,r%n
568	if(index (p%mag%name,name(1)(1:len_trim(name(1))) ) /=0) then
569	call EL_POL_force(p,poly(1))
570	dname(1)=.true.
571	elseif(index(p%mag%name,name(2)(1:len_trim(name(2))))/=0) then
572	call EL_POL_force(p,poly(2))
573	dname(2)=.true.
574	endif
575
576	p=>p%next
577	enddo
578
579	if(.not.(dname(1).and.dname(2))) then
580	CALL ELP_TO_EL(R)
581	write(6,*) " lattice_fit_TUNE_gmap_auto ---> FAILED"
582	return
583	endif
584
585
586	epsr=abs(epsf)
587
588	np=2
589
590	allocate(eq(neq))
591
592	nt=neq+np
593	STATE=((((my_state+nocavity0)-delta0)+only_4d0)-RADIATION0)
594
595	CALL INIT(STATE,no,NP)
596
597
598	! DO I=1,NPOLY
599	! R=POLY(i)
600	! ENDDO
601
602
603	CLOSED(:)=0.0_dp
604	it=0
605	100 continue
606	it=it+1
607
608	CALL FIND_ORBIT(R,CLOSED,1,STATE,1e-5_dp)
609	write(6,*) "closed orbit "
610	write(6,*) CLOSED
611
612
613	CALL INIT(STATE,no,NP)
614	CALL ALLOC(NORM)
615	CALL ALLOC(Y)
616	CALL ALLOC(EQ)
617	call alloc(id)
618
619	id=1
620	Y=CLOSED+id
621
622	CALL TRACK(R,Y,1,+STATE)
623	write(6,*) "c_%no,c_%nv,c_%nd,c_%nd2"
624	write(6,*) c_%no,c_%nv,c_%nd,c_%nd2
625	write(6,*) "c_%ndpt,c_%npara,c_%npara,c_%np_pol"
626	write(6,*) c_%ndpt,c_%npara,c_%npara,c_%np_pol
627	NORM=Y
628	gam(1)=(norm%a_t%v(2).sub.'1')2+(norm%a_t%v(2).sub.'01')2
629	gam(2)=(norm%a_t%v(4).sub.'001')2+(norm%a_t%v(4).sub.'0001')2
630	write(6,*) " Gamma= ",GAM
631	! CALL KANALNUMMER(MF)
632	OPEN(UNIT=1111,FILE='GAMMA.TXT')
633	WRITE(1111,*) " Gamma= ",GAM
634
635	write(6,*) " tunes ",NORM%TUNE(1), NORM%TUNE(2)
636
637	eq(1)= ((NORM%dhdj%v(1)).par.'0000')-targ(1)
638	eq(2)= ((NORM%dhdj%v(2)).par.'0000')-targ(2)
639	epsnow=abs(eq(1))+abs(eq(2))
640	call kanalnummer(SCRATCHFILE)
641	OPEN(UNIT=SCRATCHFILE,FILE='EQUATION.TXT')
642	rewind scratchfile
643
644	do i=1,neq
645	eq(i)=eq(i)<=c_%npara
646	enddo
647	do i=1,neq
648	call daprint(eq(i),scratchfile)
649	enddo
650	close(SCRATCHFILE)
651	CALL KILL(NORM)
652	CALL KILL(Y)
653	CALL KILL(id)
654	CALL KILL(EQ)
655
656
657
658	CALL INIT(1,nt)
659	call alloc(g,nt)
660	call kanalnummer(SCRATCHFILE)
661	OPEN(UNIT=SCRATCHFILE,FILE='EQUATION.TXT')
662	rewind scratchfile
663	do i=np+1,nt
664	call read(g%v(i),scratchfile)
665	enddo
666	close(SCRATCHFILE)
667
668	call alloc(t)
669	do i=1,np
670	g%v(i)=1.0_dp.mono.i
671	do j=np+1,nt
672	t=g%v(j).d.i
673	g%v(i)=g%v(i)+(1.0_dp.mono.j)*t
674	enddo
675	enddo
676	CALL KILL(t)
677
678	g=g.oo.(-1)
679	tpsafit=0.0_dp
680	tpsafit(1:nt)=g
681
682	SET_TPSAFIT=.true.
683
684	p=>r%start
685	do i=1,r%n
686	if(index (p%mag%name,name(1)(1:len_trim(name(1))) ) /=0) then
687	call EL_POL_force(p,poly(1))
688	elseif(index(p%mag%name,name(2)(1:len_trim(name(2))))/=0) then
689	call EL_POL_force(p,poly(2))
690	endif
691
692	p=>p%next
693	enddo
694
695	SET_TPSAFIT=.false.
696
697	CALL ELP_TO_EL(R)
698
699	! write(6,*) " more "
700	! read(5,*) more
701	if(it>=max_fit_iter) goto 101
702	if(epsnow<=epsr) goto 102
703	GOTO 100
704
705	101 continue
706	write(6,*) " warning did not converge "
707
708	102 continue
709	CALL KILL_PARA(R)
710	deallocate(eq)
711
712	end subroutine lattice_fit_TUNE_gmap_auto
713
714	subroutine lattice_fit_TUNE_gmap(R,my_state,EPSF,POLY,NPOLY,TARG,NP)
715	IMPLICIT NONE
716	TYPE(layout), target,intent(inout):: R
717	TYPE(POL_BLOCK), intent(inout),dimension(:)::POLY
718	INTEGER, intent(in):: NPOLY,NP
719	real(dp) , intent(IN),dimension(:)::TARG
720	real(dp) CLOSED(6)
721	TYPE(INTERNAL_STATE), intent(IN):: my_STATE
722	TYPE(INTERNAL_STATE) STATE
723	INTEGER I,SCRATCHFILE, MF
724	TYPE(TAYLOR), allocatable:: EQ(:)
725	TYPE(REAL_8) Y(6)
726	TYPE(NORMALFORM) NORM
727	integer :: neq=2, no=2,nt,j,it
728	type(damap) id
729	type(gmap) g
730	TYPE(TAYLOR)t
731	real(dp) epsf,epsr,epsnow,gam(2)
732	! EPSF=.0001
733	epsr=abs(epsf)
734
735	allocate(eq(neq))
736
737	nt=neq+np
738	STATE=((((my_state+nocavity0)-delta0)+only_4d0)-RADIATION0)
739
740	CALL INIT(STATE,no,NP,BERZ)
741
742	SET_TPSAFIT=.FALSE.
743
744	DO I=1,NPOLY
745	R=POLY(i)
746	ENDDO
747	CLOSED(:)=0.0_dp
748	it=0
749	100 continue
750	it=it+1
751
752	CALL FIND_ORBIT(R,CLOSED,1,STATE,1e-5_dp)
753	write(6,*) "closed orbit "
754	write(6,*) CLOSED
755
756
757	CALL INIT(STATE,no,NP,BERZ)
758	CALL ALLOC(NORM)
759	CALL ALLOC(Y)
760	CALL ALLOC(EQ)
761	call alloc(id)
762
763	id=1
764	Y=CLOSED+id
765
766	CALL TRACK(R,Y,1,+STATE)
767	write(6,*) "c_%no,c_%nv,c_%nd,c_%nd2"
768	write(6,*) c_%no,c_%nv,c_%nd,c_%nd2
769	write(6,*) "c_%ndpt,c_%npara,c_%npara,c_%np_pol"
770	write(6,*) c_%ndpt,c_%npara,c_%npara,c_%np_pol
771	NORM=Y
772	gam(1)=(norm%a_t%v(2).sub.'1')2+(norm%a_t%v(2).sub.'01')2
773	gam(2)=(norm%a_t%v(4).sub.'001')2+(norm%a_t%v(4).sub.'0001')2
774	write(6,*) " Gamma= ",GAM
775	! CALL KANALNUMMER(MF)
776	OPEN(UNIT=1111,FILE='GAMMA.TXT')
777	WRITE(1111,*) " Gamma= ",GAM
778
779	write(6,*) " tunes ",NORM%TUNE(1), NORM%TUNE(2)
780
781	eq(1)= ((NORM%dhdj%v(1)).par.'0000')-targ(1)
782	eq(2)= ((NORM%dhdj%v(2)).par.'0000')-targ(2)
783	epsnow=abs(eq(1))+abs(eq(2))
784	call kanalnummer(SCRATCHFILE)
785	OPEN(UNIT=SCRATCHFILE,FILE='EQUATION.TXT')
786	rewind scratchfile
787
788	do i=1,neq
789	eq(i)=eq(i)<=c_%npara
790	enddo
791	do i=1,neq
792	call daprint(eq(i),scratchfile)
793	enddo
794	close(SCRATCHFILE)
795	CALL KILL(NORM)
796	CALL KILL(Y)
797	CALL KILL(id)
798	CALL KILL(EQ)
799
800
801
802	CALL INIT(1,nt)
803	call alloc(g,nt)
804	call kanalnummer(SCRATCHFILE)
805	OPEN(UNIT=SCRATCHFILE,FILE='EQUATION.TXT')
806	rewind scratchfile
807	do i=np+1,nt
808	call read(g%v(i),scratchfile)
809	enddo
810	close(SCRATCHFILE)
811
812	call alloc(t)
813	do i=1,np
814	g%v(i)=1.0_dp.mono.i
815	do j=np+1,nt
816	t=g%v(j).d.i
817	g%v(i)=g%v(i)+(1.0_dp.mono.j)*t
818	enddo
819	enddo
820	CALL KILL(t)
821
822	g=g.oo.(-1)
823	tpsafit=0.0_dp
824	tpsafit(1:nt)=g
825
826	SET_TPSAFIT=.true.
827
828	DO I=1,NPOLY
829	R=POLY(i)
830	ENDDO
831	SET_TPSAFIT=.false.
832
833	CALL ELP_TO_EL(R)
834
835	! write(6,*) " more "
836	! read(5,*) more
837	if(it>=max_fit_iter) goto 101
838	if(epsnow<=epsr) goto 102
839	GOTO 100
840
841	101 continue
842	write(6,*) " warning did not converge "
843
844	102 continue
845	CALL KILL_PARA(R)
846	deallocate(eq)
847
848	end subroutine lattice_fit_TUNE_gmap
849
850	subroutine lattice_fit_CHROM_gmap(R,my_state,EPSF,POLY,NPOLY,TARG,NP)
851	IMPLICIT NONE
852	TYPE(layout),target, intent(inout):: R
853	TYPE(POL_BLOCK), intent(inout),dimension(:)::POLY
854	INTEGER, intent(in):: NPOLY,NP
855	real(dp) , intent(IN),dimension(:)::TARG
856	real(dp) CLOSED(6)
857	TYPE(INTERNAL_STATE), intent(IN):: my_STATE
858	TYPE(INTERNAL_STATE) STATE
859	INTEGER I,SCRATCHFILE
860	TYPE(TAYLOR), allocatable:: EQ(:)
861	TYPE(REAL_8) Y(6)
862	TYPE(NORMALFORM) NORM
863	integer :: neq=2, no=3,nt,j,it
864	type(damap) id
865	type(gmap) g
866	TYPE(TAYLOR)t
867	real(dp) epsf,epsr,epsnow,CHROM(2)
868	! EPSF=.0001
869	epsr=abs(epsf)
870
871	allocate(eq(neq))
872
873	nt=neq+np
874	STATE=((((my_state+nocavity0)+delta0)+only_4d0)-RADIATION0)
875
876	CALL INIT(STATE,no,NP,BERZ)
877
878	SET_TPSAFIT=.FALSE.
879
880	DO I=1,NPOLY
881	R=POLY(i)
882	ENDDO
883	CLOSED(:)=0.0_dp
884	it=0
885	100 continue
886	it=it+1
887
888	CALL FIND_ORBIT(R,CLOSED,1,STATE,1e-5_dp)
889	write(6,*) "closed orbit ", CHECK_STABLE
890	write(6,*) CLOSED
891
892
893	CALL INIT(STATE,no,NP,BERZ)
894	CALL ALLOC(NORM)
895	CALL ALLOC(Y)
896	CALL ALLOC(EQ)
897	call alloc(id)
898
899	id=1
900	Y=CLOSED+id
901
902	CALL TRACK(R,Y,1,+STATE)
903	NORM=Y
904	write(6,*) " tunes ",NORM%TUNE(1), NORM%TUNE(2), CHECK_STABLE
905	CHROM(1)=(NORM%dhdj%v(1)).SUB.'00001'
906	CHROM(2)=(NORM%dhdj%v(2)).SUB.'00001'
907	write(6,*) " CHROM ",CHROM
908
909	eq(1)= ((NORM%dhdj%v(1)).par.'00001')-targ(1)
910	eq(2)= ((NORM%dhdj%v(2)).par.'00001')-targ(2)
911	epsnow=abs(eq(1))+abs(eq(2))
912	call kanalnummer(SCRATCHFILE)
913	OPEN(UNIT=SCRATCHFILE,FILE='EQUATION.TXT')
914	rewind scratchfile
915
916	do i=1,neq
917	eq(i)=eq(i)<=c_%npara
918	enddo
919	do i=1,neq
920	call daprint(eq(i),scratchfile)
921	enddo
922	close(SCRATCHFILE)
923	CALL KILL(NORM)
924	CALL KILL(Y)
925	CALL KILL(id)
926	CALL KILL(EQ)
927
928
929
930	CALL INIT(1,nt)
931	call alloc(g,nt)
932	call kanalnummer(SCRATCHFILE)
933	OPEN(UNIT=SCRATCHFILE,FILE='EQUATION.TXT')
934	rewind scratchfile
935	do i=np+1,nt
936	call read(g%v(i),scratchfile)
937	enddo
938	close(SCRATCHFILE)
939
940	call alloc(t)
941	do i=1,np
942	g%v(i)=1.0_dp.mono.i
943	do j=np+1,nt
944	t=g%v(j).d.i
945	g%v(i)=g%v(i)+(1.0_dp.mono.j)*t
946	enddo
947	enddo
948	CALL KILL(t)
949
950	g=g.oo.(-1)
951	tpsafit=0.0_dp
952	tpsafit(1:nt)=g
953
954	SET_TPSAFIT=.true.
955
956	DO I=1,NPOLY
957	R=POLY(i)
958	ENDDO
959	SET_TPSAFIT=.false.
960
961	CALL ELP_TO_EL(R)
962
963	! write(6,*) " more "
964	! read(5,*) more
965	if(it>=max_fit_iter) goto 101
966	if(epsnow<=epsr) goto 102
967	GOTO 100
968
969	101 continue
970	write(6,*) " warning did not converge "
971
972	102 continue
973	CALL KILL_PARA(R)
974	deallocate(eq)
975
976	end subroutine lattice_fit_CHROM_gmap
977
978	subroutine lattice_fit_CHROM_gmap2(R,my_state,EPSF,POLY,NPOLY,TARG,np,n_extra,mf)
979	IMPLICIT NONE
980	integer ipause, mypause
981	TYPE(layout),target, intent(inout):: R
982	TYPE(POL_BLOCK), intent(inout),dimension(:)::POLY
983	INTEGER, intent(in):: NPOLY,np
984	real(dp) , intent(IN),dimension(:)::TARG
985	real(dp) CLOSED(6),co
986	TYPE(INTERNAL_STATE), intent(IN):: my_STATE
987	TYPE(INTERNAL_STATE) STATE
988	INTEGER I,SCRATCHFILE
989	TYPE(TAYLOR), allocatable:: EQ(:)
990	TYPE(REAL_8) Y(6)
991	TYPE(NORMALFORM) NORM
992	integer :: neq,no=3,nt,j,it,n_extra,mf
993	type(damap) id
994	type(vecresonance) vr
995	type(pbresonance) fr
996	type(gmap) g
997	TYPE(TAYLOR)t
998	real(dp) epsf,epsr,epsnow,CHROM(2)
999	integer, allocatable:: res(:,:),jt(:),je(:)
1000	real(dp), allocatable :: mat(:,:),v0(:)
1001	integer kb,kbmax,kpos,npi,ier
1002
1003
1004	neq=2
1005	allocate(res(n_extra,4))
1006	allocate(jt(5))
1007	res=0
1008	do i=1,n_extra
1009	read(mf,*) res(i,:)
1010	enddo
1011	! EPSF=.0001
1012	epsr=abs(epsf)
1013	neq=neq+2*n_extra
1014
1015	allocate(eq(neq))
1016
1017	nt=neq+np
1018	allocate(mat(nt,nt),v0(nt))
1019
1020	kbmax=0
1021
1022	DO I=1,NPOLY
1023	if(POLY(i)%nb>kbmax) kbmax= POLY(i)%nb
1024	ENDDO
1025	SET_TPSAFIT=.FALSE.
1026
1027	DO I=1,NPOLY
1028	R=POLY(i)
1029	ENDDO
1030
1031
1032	it=0
1033	100 continue
1034	it=it+1
1035	mat=0.0_dp
1036	v0=0.0_dp
1037	do i=1,np
1038	mat(i,i)=1.0_dp
1039	enddo
1040
1041	do kb=1,kbmax
1042	STATE=((((my_state+nocavity0)+delta0)+only_4d0)-RADIATION0)
1043
1044	DO I=1,NPOLY
1045	if(POLY(i)%nb==kb) then
1046	npi=POLY(i)%np
1047	kpos=POLY(i)%g-1
1048	exit
1049	endif
1050	ENDDO
1051
1052	write(6,*) " np in batch ",kb," = ",npi
1053	! pause 1
1054
1055	CLOSED(:)=0.0_dp
1056
1057	CALL FIND_ORBIT(R,CLOSED,1,STATE,1e-5_dp)
1058	write(6,*) "closed orbit ", CHECK_STABLE
1059	write(6,*) CLOSED
1060
1061
1062	CALL INIT(STATE,no,NPi,BERZ)
1063	nb_=kb
1064	CALL ALLOC(NORM)
1065	CALL ALLOC(Y)
1066	CALL ALLOC(EQ)
1067	call alloc(id)
1068	call alloc(vr)
1069	call alloc(fr)
1070
1071	id=1
1072	Y=CLOSED+id
1073
1074	CALL TRACK(R,Y,1,+STATE)
1075	NORM=Y
1076	vr=norm%a%nonlinear
1077	fr=norm%a%pb
1078	! call print(vr%cos%v(2),6)
1079	! call print(vr%sin%v(2),6)
1080	! pause 1
1081	! call print(fr%cos,6)
1082	! call print(fr%sin,6)
1083	! pause 2
1084
1085
1086	write(6,*) " tunes ",NORM%TUNE(1), NORM%TUNE(2), CHECK_STABLE
1087	CHROM(1)=(NORM%dhdj%v(1)).SUB.'00001'
1088	CHROM(2)=(NORM%dhdj%v(2)).SUB.'00001'
1089	write(6,*) " CHROM ",CHROM
1090
1091	eq(1)= ((NORM%dhdj%v(1)).par.'00001')-targ(1)
1092	eq(2)= ((NORM%dhdj%v(2)).par.'00001')-targ(2)
1093	do i=1,n_extra
1094	jt=0
1095	jt(1:4)=res(i,:)
1096	jt(1)=jt(1)-1
1097	eq(2+2*i-1)= ((vr%cos%v(2)).par.jt)
1098	eq(2+2*i)= ((vr%sin%v(2)).par.jt)
1099	enddo
1100
1101	epsnow=0.0_dp
1102	do i=1,neq
1103	epsnow=abs(eq(i))+epsnow
1104	if(kb==1) then
1105	co=eq(i)
1106	if(i<=2) then
1107	co=co+targ(i)
1108	endif
1109	write(6,*) i,co
1110	endif
1111	enddo
1112	write(6,*) "epsnow ", epsnow
1113	ipause=mypause(123)
1114	call kanalnummer(SCRATCHFILE)
1115	OPEN(UNIT=SCRATCHFILE,FILE='EQUATION.TXT')
1116	rewind scratchfile
1117	allocate(je(c_%nv))
1118	! write(6,*) nt,neq,np,kpos,npi
1119	je=0
1120	do i=1,neq
1121	eq(i)=eq(i)<=c_%npara
1122	v0(i+np)=-(eq(i).sub.'0')
1123	do j=1,npi
1124	je(j)=1
1125	co=eq(i).sub.je
1126	mat(np+i,j+kpos)=co
1127	mat(j+kpos,np+i)=co
1128	je(j)=0
1129	enddo
1130	enddo
1131	deallocate(je)
1132	do i=1,neq
1133	call daprint(eq(i),scratchfile)
1134	enddo
1135	close(SCRATCHFILE)
1136	CALL KILL(NORM)
1137	CALL KILL(Y)
1138	CALL KILL(id)
1139	CALL KILL(vr)
1140	CALL KILL(fr)
1141	CALL KILL(EQ)
1142	! pause 888
1143	enddo ! kbmax
1144	write(6,*) "Iteration # ",it
1145
1146	! pause 2000
1147	! do i=1,nt
1148	! do j=1,nt
1149	! if(mat(i,j)/=zero) write(6,*) i,j,mat(i,j)
1150	! enddo
1151	! enddo
1152	call matinv(mat,mat,nt,nt,ier)
1153	if(ier/=0 ) then
1154	write(6,*) ier
1155	write(6,*) " inversion error "
1156
1157	stop
1158	endif
1159	! pause 2001
1160
1161	v0=matmul(mat,v0)
1162	tpsafit=0.0_dp
1163	tpsafit(1:nt)=v0
1164
1165	SET_TPSAFIT=.true.
1166
1167	DO I=1,NPOLY
1168	R=POLY(i)
1169	ENDDO
1170	SET_TPSAFIT=.false.
1171
1172	CALL ELP_TO_EL(R)
1173
1174	! write(6,*) " more "
1175	! read(5,*) more
1176	if(it>=max_fit_iter) goto 101
1177	if(epsnow<=epsr) goto 102
1178	GOTO 100
1179
1180	101 continue
1181	write(6,*) " warning did not converge "
1182
1183	102 continue
1184	CALL KILL_PARA(R)
1185	deallocate(mat)
1186	deallocate(eq)
1187	deallocate(res)
1188	deallocate(jt)
1189
1190	end subroutine lattice_fit_CHROM_gmap2
1191
1192	subroutine lattice_fit_CHROM_gmap1(R,my_state,EPSF,POLY,NPOLY,TARG,NP,n_extra,mf)
1193	IMPLICIT NONE
1194	integer ipause, mypause
1195	TYPE(layout),target, intent(inout):: R
1196	TYPE(POL_BLOCK), intent(inout),dimension(:)::POLY
1197	INTEGER, intent(in):: NPOLY,NP
1198	real(dp) , intent(IN),dimension(:)::TARG
1199	real(dp) CLOSED(6)
1200	TYPE(INTERNAL_STATE), intent(IN):: my_STATE
1201	TYPE(INTERNAL_STATE) STATE
1202	INTEGER I,SCRATCHFILE
1203	TYPE(TAYLOR), allocatable:: EQ(:)
1204	TYPE(REAL_8) Y(6)
1205	TYPE(NORMALFORM) NORM
1206	integer :: neq=2, no=3,nt,j,it,n_extra,mf
1207	type(damap) id
1208	type(vecresonance) vr
1209	type(pbresonance) fr
1210	type(gmap) g
1211	TYPE(TAYLOR)t
1212	real(dp) epsf,epsr,epsnow,CHROM(2)
1213	integer, allocatable:: res(:,:),jt(:)
1214
1215	neq=2
1216
1217	allocate(res(n_extra,4))
1218	allocate(jt(5))
1219	res=0
1220	do i=1,n_extra
1221	read(mf,*) res(i,:)
1222	enddo
1223	! EPSF=.0001
1224	epsr=abs(epsf)
1225	neq=neq+2*n_extra
1226	allocate(eq(neq))
1227
1228	nt=neq+np
1229	STATE=((((my_state+nocavity0)+delta0)+only_4d0)-RADIATION0)
1230
1231	! CALL INIT(STATE,no,NP,BERZ)
1232
1233	! SET_TPSAFIT=.FALSE.
1234
1235	DO I=1,NPOLY
1236	R=POLY(i)
1237	ENDDO
1238	CLOSED(:)=0.0_dp
1239	it=0
1240	100 continue
1241	it=it+1
1242
1243	CALL FIND_ORBIT(R,CLOSED,1,STATE,1e-5_dp)
1244	write(6,*) "closed orbit ", CHECK_STABLE
1245	write(6,*) CLOSED
1246
1247
1248	CALL INIT(STATE,no,NP,BERZ)
1249	CALL ALLOC(NORM)
1250	CALL ALLOC(Y)
1251	CALL ALLOC(EQ)
1252	call alloc(id)
1253	call alloc(vr)
1254	call alloc(fr)
1255
1256	id=1
1257	Y=CLOSED+id
1258
1259	CALL TRACK(R,Y,1,+STATE)
1260	NORM=Y
1261	vr=norm%a%nonlinear
1262	fr=norm%a%pb
1263	! call print(vr%cos%v(2),6)
1264	! call print(vr%sin%v(2),6)
1265	! pause 1
1266	call print(fr%cos,6)
1267	call print(fr%sin,6)
1268	ipause=mypause(2)
1269
1270
1271	write(6,*) " tunes ",NORM%TUNE(1), NORM%TUNE(2), CHECK_STABLE
1272	CHROM(1)=(NORM%dhdj%v(1)).SUB.'00001'
1273	CHROM(2)=(NORM%dhdj%v(2)).SUB.'00001'
1274	write(6,*) " CHROM ",CHROM
1275
1276	eq(1)= ((NORM%dhdj%v(1)).par.'00001')-targ(1)
1277	eq(2)= ((NORM%dhdj%v(2)).par.'00001')-targ(2)
1278	do i=1,n_extra
1279	jt=0
1280	jt(1:4)=res(i,:)
1281	jt(1)=jt(1)-1
1282	eq(2+2*i-1)= ((vr%cos%v(2)).par.jt)
1283	eq(2+2*i)= ((vr%sin%v(2)).par.jt)
1284	enddo
1285
1286	epsnow=0.0_dp
1287	do i=1,neq
1288	epsnow=abs(eq(i))+epsnow
1289	enddo
1290	call kanalnummer(SCRATCHFILE)
1291	OPEN(UNIT=SCRATCHFILE,FILE='EQUATION.TXT')
1292	rewind scratchfile
1293
1294	do i=1,neq
1295	eq(i)=eq(i)<=c_%npara
1296	enddo
1297	do i=1,neq
1298	call daprint(eq(i),scratchfile)
1299	enddo
1300	close(SCRATCHFILE)
1301	CALL KILL(NORM)
1302	CALL KILL(Y)
1303	CALL KILL(id)
1304	CALL KILL(vr)
1305	CALL KILL(fr)
1306	CALL KILL(EQ)
1307
1308
1309
1310	CALL INIT(1,nt)
1311	call alloc(g,nt)
1312	call kanalnummer(SCRATCHFILE)
1313	OPEN(UNIT=SCRATCHFILE,FILE='EQUATION.TXT')
1314	rewind scratchfile
1315	do i=np+1,nt
1316	call read(g%v(i),scratchfile)
1317	enddo
1318	close(SCRATCHFILE)
1319
1320	call alloc(t)
1321	do i=1,np
1322	g%v(i)=1.0_dp.mono.i
1323	do j=np+1,nt
1324	t=g%v(j).d.i
1325	g%v(i)=g%v(i)+(1.0_dp.mono.j)*t
1326	enddo
1327	enddo
1328	CALL KILL(t)
1329
1330	g=g.oo.(-1)
1331	tpsafit(1:nt)=g
1332
1333	SET_TPSAFIT=.true.
1334
1335	DO I=1,NPOLY
1336	R=POLY(i)
1337	ENDDO
1338	SET_TPSAFIT=.false.
1339
1340	CALL ELP_TO_EL(R)
1341
1342	! write(6,*) " more "
1343	! read(5,*) more
1344	ipause=mypause(777)
1345	if(it>=max_fit_iter) goto 101
1346	if(epsnow<=epsr) goto 102
1347	GOTO 100
1348
1349	101 continue
1350	write(6,*) " warning did not converge "
1351
1352	102 continue
1353	CALL KILL_PARA(R)
1354	deallocate(eq)
1355	deallocate(res)
1356	deallocate(jt)
1357
1358	end subroutine lattice_fit_CHROM_gmap1
1359
1360	subroutine lattice_fit_tune_CHROM_gmap(R,my_state,EPSF,POLY,NPOLY,TARG,NP)
1361	IMPLICIT NONE
1362	TYPE(layout),target, intent(inout):: R
1363	TYPE(POL_BLOCK), intent(inout),dimension(:)::POLY
1364	INTEGER, intent(in):: NPOLY,NP
1365	real(dp) , intent(IN),dimension(:)::TARG
1366	real(dp) CLOSED(6)
1367	TYPE(INTERNAL_STATE), intent(IN):: my_STATE
1368	TYPE(INTERNAL_STATE) STATE
1369	INTEGER I,SCRATCHFILE, MF
1370	TYPE(TAYLOR), allocatable:: EQ(:)
1371	TYPE(REAL_8) Y(6)
1372	TYPE(NORMALFORM) NORM
1373	integer :: neq=4, no=3,nt,j,it
1374	type(damap) id
1375	type(gmap) g
1376	TYPE(TAYLOR)t
1377	real(dp) epsf,epsr,epsnow,tune(2),CHROM(2),gam(2)
1378	! EPSF=.0001
1379	epsr=abs(epsf)
1380
1381	allocate(eq(neq))
1382
1383	nt=neq+np
1384	STATE=((((my_state+nocavity0)+delta0)+only_4d0)-RADIATION0)
1385
1386	CALL INIT(STATE,no,NP,BERZ)
1387
1388	SET_TPSAFIT=.FALSE.
1389
1390	DO I=1,NPOLY
1391	R=POLY(i)
1392	ENDDO
1393	CLOSED(:)=0.0_dp
1394	it=0
1395	100 continue
1396	it=it+1
1397
1398	CALL FIND_ORBIT(R,CLOSED,1,STATE,1e-5_dp)
1399	write(6,*) "closed orbit ", CHECK_STABLE
1400	write(6,*) CLOSED
1401
1402
1403	CALL INIT(STATE,no,NP,BERZ)
1404	CALL ALLOC(NORM)
1405	CALL ALLOC(Y)
1406	CALL ALLOC(EQ)
1407	call alloc(id)
1408
1409	id=1
1410	Y=CLOSED+id
1411
1412	CALL TRACK(R,Y,1,+STATE)
1413	NORM=Y
1414	gam(1)=(norm%a_t%v(2).sub.'1')2+(norm%a_t%v(2).sub.'01')2
1415	gam(2)=(norm%a_t%v(4).sub.'001')2+(norm%a_t%v(4).sub.'0001')2
1416	write(6,*) " Gamma= ",GAM
1417	! CALL KANALNUMMER(MF)
1418	OPEN(UNIT=1111,FILE='GAMMA.TXT')
1419	WRITE(1111,*) " Gamma= ",GAM
1420
1421	write(6,*) " tunes ",NORM%TUNE(1), NORM%TUNE(2), CHECK_STABLE
1422	tune(1)=(NORM%dhdj%v(1)).SUB.'0000'
1423	tune(2)=(NORM%dhdj%v(2)).SUB.'0000'
1424	CHROM(1)=(NORM%dhdj%v(1)).SUB.'00001'
1425	CHROM(2)=(NORM%dhdj%v(2)).SUB.'00001'
1426	write(6,*) " CHROM ",CHROM
1427
1428	eq(1)= ((NORM%dhdj%v(1)).par.'00000')-targ(1)
1429	eq(2)= ((NORM%dhdj%v(2)).par.'00000')-targ(2)
1430	eq(3)= ((NORM%dhdj%v(1)).par.'00001')-targ(3)
1431	eq(4)= ((NORM%dhdj%v(2)).par.'00001')-targ(4)
1432	epsnow=abs(eq(1))+abs(eq(2))+abs(eq(3))+abs(eq(4))
1433	call kanalnummer(SCRATCHFILE)
1434	OPEN(UNIT=SCRATCHFILE,FILE='EQUATION.TXT')
1435	rewind scratchfile
1436
1437	do i=1,neq
1438	eq(i)=eq(i)<=c_%npara
1439	enddo
1440	do i=1,neq
1441	call daprint(eq(i),scratchfile)
1442	enddo
1443	close(SCRATCHFILE)
1444	CALL KILL(NORM)
1445	CALL KILL(Y)
1446	CALL KILL(id)
1447	CALL KILL(EQ)
1448
1449
1450
1451	CALL INIT(1,nt)
1452	call alloc(g,nt)
1453	call kanalnummer(SCRATCHFILE)
1454	OPEN(UNIT=SCRATCHFILE,FILE='EQUATION.TXT')
1455	rewind scratchfile
1456	do i=np+1,nt
1457	call read(g%v(i),scratchfile)
1458	enddo
1459	close(SCRATCHFILE)
1460
1461	call alloc(t)
1462	do i=1,np
1463	g%v(i)=1.0_dp.mono.i
1464	do j=np+1,nt
1465	t=g%v(j).d.i
1466	g%v(i)=g%v(i)+(1.0_dp.mono.j)*t
1467	enddo
1468	enddo
1469	CALL KILL(t)
1470
1471	g=g.oo.(-1)
1472	tpsafit(1:nt)=g
1473
1474	SET_TPSAFIT=.true.
1475
1476	DO I=1,NPOLY
1477	R=POLY(i)
1478	ENDDO
1479	SET_TPSAFIT=.false.
1480
1481	CALL ELP_TO_EL(R)
1482
1483	! write(6,*) " more "
1484	! read(5,*) more
1485	if(it>=max_fit_iter) goto 101
1486	if(epsnow<=epsr) goto 102
1487	GOTO 100
1488
1489	101 continue
1490	write(6,*) " warning did not converge "
1491
1492	102 continue
1493	CALL KILL_PARA(R)
1494	deallocate(eq)
1495
1496	end subroutine lattice_fit_tune_CHROM_gmap
1497
1498
1499	subroutine lattice_PRINT_RES_FROM_A(R,my_state,NO,EMIT0,MRES,FILENAME)
1500	IMPLICIT NONE
1501	TYPE(layout),target, intent(inout):: R
1502	real(dp) CLOSED(6)
1503	TYPE(INTERNAL_STATE), intent(IN):: my_STATE
1504	TYPE(INTERNAL_STATE) STATE
1505	TYPE(REAL_8) Y(6)
1506	TYPE(NORMALFORM) NORM
1507	type(damap) id
1508	TYPE(PBresonance)Hr
1509	CHARACTER(*) FILENAME
1510	REAL(DP) PREC,EMIT0(2),STR
1511	INTEGER NO,MF,MRES(4)
1512
1513	PREC=1e-10_dp
1514
1515	STATE=((((my_state+nocavity0)-delta0)+only_4d0)-RADIATION0)
1516
1517
1518
1519	CALL FIND_ORBIT(R,CLOSED,1,STATE,1e-5_dp)
1520	write(6,*) "closed orbit "
1521	write(6,*) CLOSED
1522
1523
1524	CALL INIT(STATE,no,0,BERZ)
1525
1526	CALL ALLOC(NORM)
1527	CALL ALLOC(Y)
1528	CALL ALLOC(hr)
1529	call alloc(id) ! LOOK AT PAGE 143 OF THE BOOK
1530	id=1
1531	Y=CLOSED+id
1532
1533	CALL TRACK(R,Y,1,STATE)
1534
1535
1536	NORM=Y
1537
1538
1539	hr=NORM%a%pb
1540
1541	CALL KANALNUMMER(MF)
1542
1543	OPEN(UNIT=MF,FILE=FILENAME)
1544
1545	CALL PRINT(HR%COS,6,PREC)
1546	CALL PRINT(HR%SIN,6,PREC)
1547	STR=(HR%COS%H.SUB.MRES)2+(HR%SIN%H.SUB.MRES)2; STR=SQRT(STR)
1548	WRITE(6,*) "RESONANCE = ",MRES,STR
1549	WRITE(MF,*) "RESONANCE = ",MRES,STR
1550	CALL PRINT(HR%COS,MF,PREC)
1551	CALL PRINT(HR%SIN,MF,PREC)
1552	WRITE(MF,*) " SCALE AT EMIT = ",EMIT0
1553	ID=1
1554	ID%V(1)=ID%V(1)*SQRT(EMIT0(1))
1555	ID%V(2)=ID%V(2)*SQRT(EMIT0(1))
1556	ID%V(3)=ID%V(3)*SQRT(EMIT0(2))
1557	ID%V(4)=ID%V(4)*SQRT(EMIT0(2))
1558	HR%COS%H=HR%COS%H*ID
1559	HR%SIN%H=HR%SIN%H*ID
1560	CALL PRINT(HR%COS,MF,PREC)
1561	CALL PRINT(HR%SIN,MF,PREC)
1562
1563
1564	CLOSE(MF)
1565	CALL KILL(NORM)
1566	CALL KILL(Y)
1567	CALL KILL(hr)
1568
1569	end subroutine lattice_PRINT_RES_FROM_A
1570
1571	! call lattice_random_error(my_ring,name,i1,i2,cut,n,addi,integrated,cn,cns,sc)
1572
1573
1574
1575	subroutine lattice_random_error_new(R,nom,full,iseed,cut,n,addi,integrated,cn,cns,per,pr)
1576	use gauss_dis
1577	IMPLICIT NONE
1578	TYPE(layout),target, intent(inout):: R
1579	integer n,addi,ic,i,iseed,j
1580	character(nlp) nom
1581	type(fibre), pointer :: p
1582	logical(lp) integrated,f1,f2,full,pr
1583	real(dp) x,bn,cn,cns,cut,per
1584
1585	if(iseed/=0) call gaussian_seed(iseed)
1586
1587
1588	call context(nom)
1589
1590	ic=0
1591	p=>r%start
1592	do i=1,r%n
1593	f1=.false.
1594	f2=.false.
1595	if(full) then
1596	f1=(p%mag%name ==nom )
1597	else
1598	f1=( index(p%mag%name,nom(1:len_trim(nom))) >0)
1599	endif
1600
1601
1602	if(f1.and.per<=0) then
1603	call GRNF(X,cut)
1604	bn=cns+cn*x
1605	if(integrated.and.p%mag%p%ld/=0.0_dp) then
1606	bn=bn/p%mag%l
1607	endif
1608	if(bn/=0.0_dp) call add(p,n,addi,bn)
1609	f2=.true.
1610	endif
1611
1612	if(f1.and.per/=0.0_dp) then
1613	call GRNF(X,cut)
1614	do j=p%mag%p%nmul,1,-1
1615	! if(n>0) then
1616	bn=p%mag%bn(j)
1617	bn=bn(1.0_dp+xper)
1618	call add(p,j,0,bn)
1619	! else
1620	bn=p%mag%an(j)
1621	bn=bn(1.0_dp+xper)
1622	! endif
1623	call add(p,-j,0,bn)
1624
1625	enddo
1626	f2=.true.
1627	endif
1628	if(f2) then
1629	ic=ic+1
1630	if(pr) write(6,*) p%mag%name
1631	endif
1632	p=>P%next
1633	enddo
1634
1635	if(pr) write(6,*) ic," Magnets modified "
1636
1637	end subroutine lattice_random_error_new
1638
1639
1640	subroutine toggle_verbose
1641	implicit none
1642	verbose=.not.verbose
1643	end subroutine toggle_verbose
1644
1645
1646
1647
1648
1649	! linked
1650
1651	SUBROUTINE FIND_ORBIT_LAYOUT(RING,FIX,LOC,STATE,TURNS) ! Finds orbit with TPSA in State or compatible state
1652	IMPLICIT NONE
1653	TYPE(layout),target,INTENT(INOUT):: RING
1654	INTEGER, OPTIONAL:: TURNS
1655	real(dp) FIX(6),DIX(6),xdix,xdix0,tiny,freq
1656	TYPE(REAL_8) X(6)
1657	TYPE(DAMAP) MX,SX,SXI,IS
1658	integer NO1,ND2,I,IU,LOC,ITE,npara
1659	TYPE(INTERNAL_STATE),optional, intent(in) :: STATE
1660	TYPE(INTERNAL_STATE) stat
1661	TYPE (fibre), POINTER :: C
1662	logical(lp) c_da,s_da
1663	INTEGER TURNS0
1664	s_da=c_%stable_da
1665	c_da=c_%check_da
1666	! APERTURE=c_%APERTURE_FLAG
1667
1668	! c_%APERTURE_FLAG=.false.
1669	c_%stable_da=.true.
1670	c_%check_da=.true.
1671	TURNS0=1
1672	freq=0.0_dp
1673	IF(PRESENT(TURNS)) TURNS0=TURNS
1674	Nullify(C);
1675	if(.not.ring%closed) then
1676	w_p=0
1677	w_p%nc=1
1678	w_p%fc='((1X,a72))'
1679	w_p%c(1)= " This line is not ring : FIND_ORBIT_LAYOUT "
1680	! call !write_e(100)
1681	endif
1682	dix(:)=0.0_dp
1683	tiny=1e-40_dp
1684	xdix0=1e4_dp*DEPS_tracking
1685	NO1=1
1686
1687	if(.not.present(STATE)) then
1688	IF(default%NOCAVITY) THEN
1689	! ND1=2
1690	stat=default+only_4d-spin0
1691	ELSE
1692	! ND1=3
1693	STAT=default-spin0
1694	C=>RING%START
1695	do i=1,RING%n
1696	if(C%magp%kind==kind4.OR.C%magp%kind==kind21) goto 101
1697	C=>C%NEXT
1698	enddo
1699	w_p=0
1700	w_p%nc=2
1701	w_p%fc='((1X,a72,/),(1X,a72))'
1702	w_p%c(1)= " No Cavity in the Line "
1703	w_p%c(2)= " FIND_ORBIT_LAYOUT will crash "
1704	! call !write_e(111)
1705	ENDIF
1706	else
1707	IF(STATE%NOCAVITY) THEN
1708	! ND1=2
1709	STAT=STATE+only_4d0-spin0
1710	ELSE
1711	! ND1=3
1712	STAT=STATE-spin0
1713	C=>RING%START
1714	do i=1,RING%n
1715	if(C%magp%kind==kind4.OR.C%magp%kind==kind21) goto 101
1716	C=>C%NEXT
1717	enddo
1718	w_p=0
1719	w_p%nc=2
1720	w_p%fc='((1X,a72,/),(1X,a72))'
1721	w_p%c(1)= " No Cavity in the Line "
1722	w_p%c(2)= " FIND_ORBIT_LAYOUT will crash "
1723	! call !write_e(112)
1724	ENDIF
1725	endif
1726	101 continue
1727	! ND2=2*ND1
1728	if(stat%totalpath==1.and.(.not.stat%nocavity)) then
1729	C=>RING%START
1730	freq=0.0_dp
1731	i=1
1732	xdix=0.0_dp
1733	do while(i<=RING%n)
1734	if(associated(c%magp%freq)) then
1735	IF(FREQ==0.0_dp) THEN
1736	freq=c%magp%freq
1737	ELSEIF(c%magp%freq/=0.0_dp.AND.c%magp%freq<FREQ) THEN
1738	freq=c%magp%freq
1739	ENDIF
1740	endif
1741	IF(stat%TIME) THEN
1742	XDIX=XDIX+c%mag%P%LD/c%BETA0
1743	ELSE
1744	XDIX=XDIX+c%mag%P%LD
1745	ENDIF
1746	c=>c%next
1747	i=i+1
1748	enddo
1749	if(freq==0.0_dp) then
1750	w_p=0
1751	w_p%nc=2
1752	w_p%fc='((1X,a72,/),(1X,a72))'
1753	w_p%c(1)= " No Cavity in the Line or Frequency = 0 "
1754	w_p%c(2)= " FIND_ORBIT_LAYOUT will crash "
1755	! call !write_e(113)
1756	endif
1757	IF(RING%HARMONIC_NUMBER>0) THEN
1758	FREQ=RING%HARMONIC_NUMBER*CLIGHT/FREQ
1759	stop 475
1760	ELSE
1761	XDIX=XDIX*FREQ/CLIGHT
1762	FREQ=NINT(XDIX)*CLIGHT/FREQ
1763	ENDIF
1764	endif
1765	CALL INIT(STAT,NO1,0,BERZ,ND2,NPARA)
1766
1767	! call init(NO1,ND1,0,0,.TRUE.)
1768
1769
1770	CALL ALLOC(X,6)
1771	CALL ALLOC(MX)
1772	CALL ALLOC(SX)
1773	CALL ALLOC(SXI)
1774	CALL ALLOC(IS)
1775
1776
1777	3 continue
1778	X=NPARA
1779	DO I=1,6
1780	X(I)=FIX(I)
1781	ENDDO
1782
1783	DO I=1,TURNS0
1784	CALL TRACK(RING,X,LOC,STAT)
1785	if(.not.check_stable.or.(.not.c_%stable_da)) then
1786	CALL KILL(X,6)
1787	CALL KILL(MX)
1788	CALL KILL(SX)
1789	CALL KILL(SXI)
1790	CALL KILL(IS)
1791	w_p=0
1792	w_p%nc=1
1793	w_p%fc='((1X,a72))'
1794	write(w_p%c(1),'(a30,i4)') " Lost in Fixed Point Searcher ",1
1795	messagelost(len_trim(messagelost)+1:255)=" -> Lost in Fixed Point Searcher "
1796	! call ! WRITE_I
1797
1798	return
1799	endif
1800	ENDDO
1801	x(6)=x(6)-TURNS0*freq
1802
1803	IS=1
1804	MX=X
1805	SX=MX-IS
1806	DIX=SX
1807	DO I=1,ND2
1808	DIX(I)=FIX(I)-DIX(I)
1809	enddo
1810	IS=0
1811	IS=DIX
1812
1813	SXI=SX**(-1)
1814	SX=SXI.o.IS
1815	DIX=SX
1816	DO I=1,ND2
1817	FIX(I)=FIX(I)+DIX(I)
1818	ENDDO
1819
1820	xdix=0.0_dp
1821	do iu=1,ND2
1822	xdix=abs(dix(iu))+xdix
1823	enddo
1824
1825	if(verbose) write(w_p%c(1),'(a22,g21.14)') " Convergence Factor = ",xdix
1826	! if(verbose) ! call ! WRITE_I
1827	if(xdix.gt.deps_tracking) then
1828	ite=1
1829	else
1830	if(xdix.ge.xdix0.or.xdix<=tiny) then
1831	ite=0
1832	else
1833	ite=1
1834	xdix0=xdix
1835	endif
1836	endif
1837	if(ite.eq.1) then
1838	GOTO 3
1839
1840	endif
1841	CALL KILL(X,6)
1842	CALL KILL(MX)
1843	CALL KILL(SX)
1844	CALL KILL(SXI)
1845	CALL KILL(IS)
1846	c_%check_da=c_da
1847	! c_%APERTURE_FLAG=APERTURE
1848	c_%stable_da=s_da
1849
1850	! write(6,*) " 2 ",APERTURE,APERTURE_FLAG
1851
1852	END SUBROUTINE FIND_ORBIT_LAYOUT
1853
1854
1855	SUBROUTINE FIND_ORBIT_LAYOUT_noda(RING,FIX,LOC,STATE,eps,TURNS) ! Finds orbit without TPSA in State or compatible state
1856	IMPLICIT NONE
1857	TYPE(layout),target,INTENT(INOUT):: RING
1858	real(dp) , intent(inOUT) :: FIX(6)
1859	INTEGER , intent(in) :: LOC
1860	INTEGER, OPTIONAL::TURNS
1861	real(dp) , optional,intent(in) :: eps
1862	TYPE(INTERNAL_STATE),optional, intent(in) :: STATE
1863	TYPE(INTERNAL_STATE) stat
1864
1865	real(dp) DIX(6),xdix,xdix0,tiny,freq
1866	real(dp) X(6),Y(6),MX(6,6),sxi(6,6),SX(6,6)
1867	integer NO1,ND2,I,IU,ITE,ier,j,ITEM
1868	TYPE (fibre), POINTER :: C
1869	logical(lp) APERTURE
1870	INTEGER TURNS0,trackflag
1871	TURNS0=1
1872	trackflag=0
1873	IF(PRESENT(TURNS)) TURNS0=TURNS
1874	freq=0.0_dp
1875	APERTURE=c_%APERTURE_FLAG
1876	c_%APERTURE_FLAG=.false.
1877	messagelost=' Orbit most likely found'
1878	if(.not.present(eps)) then
1879	if(.not.present(STATE)) then
1880	call FIND_ORBIT_LAYOUT(RING,FIX,LOC,TURNS=TURNS0)
1881	else
1882	call FIND_ORBIT_LAYOUT(RING,FIX,LOC,STATE,TURNS=TURNS0)
1883	endif
1884	c_%APERTURE_FLAG=APERTURE
1885	return
1886	endif
1887
1888
1889	Nullify(C);
1890
1891	if(.not.ring%closed) then
1892	write(6,*) " This line is not ring : FIND_ORBIT_LAYOUT_noda "
1893	check_stable=.false.
1894	! call !write_e(100)
1895	endif
1896	dix(:)=0.0_dp
1897	tiny=1e-40_dp
1898	xdix0=1e4_dp*DEPS_tracking
1899	NO1=1
1900	if(.not.present(STATE)) then
1901	IF(default%NOCAVITY) THEN
1902	! ND1=2
1903	stat=default+ only_4d-spin0
1904	ELSE
1905	! ND1=3
1906	STAT=default-spin0
1907	C=>RING%START
1908	do i=1,RING%n
1909	if(C%magp%kind==kind4.OR.C%magp%kind==kind21) goto 101
1910	C=>C%NEXT
1911	enddo
1912	messagelost= " FIND_ORBIT_LAYOUT will crash : exiting"
1913	check_stable=.false.
1914	return
1915	ENDIF
1916	else
1917	IF(STATE%NOCAVITY) THEN
1918	ND2=4
1919	STAT=STATE+only_4d0-spin0
1920	ELSE
1921	ND2=6
1922	STAT=STATE-spin0
1923	C=>RING%START
1924	do i=1,RING%n
1925	if(C%magp%kind==kind4.OR.C%magp%kind==kind21) goto 101
1926	C=>C%NEXT
1927	enddo
1928	messagelost=" State present; no cavity: FIND_ORBIT_LAYOUT will crash => exiting"
1929	check_stable=.false.
1930	return
1931	ENDIF
1932	endif
1933	101 continue
1934
1935
1936	if(stat%totalpath==1.and.(.not.stat%nocavity)) then
1937	C=>RING%START
1938	freq=0.0_dp
1939	i=1
1940	xdix=0.0_dp
1941	do while(i<=RING%n)
1942	if(associated(c%magp%freq)) then
1943	IF(FREQ==0.0_dp) THEN
1944	freq=c%magp%freq
1945	ELSEIF(c%magp%freq/=0.0_dp.AND.c%magp%freq<FREQ) THEN
1946	freq=c%magp%freq
1947	ENDIF
1948	endif
1949	XDIX=XDIX+c%mag%P%LD/c%BETA0
1950	c=>c%next
1951	i=i+1
1952	enddo
1953	if(freq==0.0_dp) then
1954
1955	messagelost= " No Cavity in the Line or Frequency = 0 (totalpath==1)"
1956	check_stable=.false.
1957	return
1958	endif
1959	IF(RING%HARMONIC_NUMBER>0) THEN
1960	FREQ=RING%HARMONIC_NUMBER*CLIGHT/FREQ
1961	ELSE
1962	XDIX=XDIX*FREQ/CLIGHT
1963	FREQ=NINT(XDIX)*CLIGHT/FREQ
1964	ENDIF
1965	endif
1966
1967
1968
1969	ITEM=0
1970	3 continue
1971	ITEM=ITEM+1
1972	X=FIX
1973
1974	DO I=1,TURNS0
1975	! CALL TRACK(RING,X,LOC,STAT)
1976	call TRACK(RING,X,LOC,STAT)
1977	if(.not.check_stable) then
1978	messagelost(len_trim(messagelost)+1:255)=" -> Unstable tracking: closed orbit not found"
1979	c_%APERTURE_FLAG=APERTURE
1980	return
1981	endif
1982	! if(.not.check_stable) then
1983	! w_p=0
1984	! w_p%nc=1
1985	! w_p%fc='((1X,a72))'
1986	! write(w_p%c(1),'(a30,i4)') " Lost in Fixed Point Searcher ",2
1987	! ! call ! WRITE_I
1988
1989	! return
1990	! endif
1991
1992	ENDDO
1993	! write(6,*) x
1994	x(6)=x(6)-freq*turns0
1995
1996	mx=0.0_dp
1997	DO J=1,ND2
1998
1999	Y=FIX
2000	Y(J)=FIX(J)+EPS
2001	DO I=1,TURNS0
2002	CALL TRACK(RING,Y,LOC,STAT)
2003	if(.not.check_stable) then
2004	messagelost(len_trim(messagelost)+1:255)=" -> Unstable while tracking small rays around the guessed orbit "
2005	! fixed_found=my_false
2006	c_%APERTURE_FLAG=APERTURE
2007	return
2008	endif
2009	ENDDO
2010	y(6)=y(6)-freq*turns0
2011
2012	do i=1,ND2
2013	MX(I,J)=(Y(i)-X(i))/eps
2014	enddo
2015
2016	ENDDO
2017
2018
2019	SX=MX;
2020	DO I=1,nd2 ! 6 before
2021	SX(I,I)=MX(I,I)-1.0_dp
2022	ENDDO
2023
2024	DO I=1,ND2
2025	DIX(I)=FIX(I)-X(I)
2026	enddo
2027
2028	CALL matinv(SX,SXI,ND2,6,ier)
2029	IF(IER==132) then
2030	messagelost= " Inversion failed in FIND_ORBIT_LAYOUT_noda"
2031	check_stable=.false.
2032	return
2033	endif
2034
2035	x=0.0_dp
2036	do i=1,nd2
2037	do j=1,nd2
2038	x(i)=sxi(i,j)*dix(j)+x(i)
2039	enddo
2040	enddo
2041	dix=x
2042	DO I=1,ND2
2043	FIX(I)=FIX(I)+DIX(I)
2044	ENDDO
2045
2046	xdix=0.0_dp
2047	do iu=1,ND2
2048	xdix=abs(dix(iu))+xdix
2049	enddo
2050	! write(6,*) " Convergence Factor = ",nd2,xdix,deps_tracking
2051	! pause 123321
2052
2053	if(verbose) write(6,*) " Convergence Factor = ",xdix
2054	! if(verbose) ! call ! WRITE_I
2055	if(xdix.gt.deps_tracking) then
2056	ite=1
2057	else
2058	if(xdix.ge.xdix0.or.xdix<=tiny) then
2059	ite=0
2060	else
2061	ite=1
2062	xdix0=xdix
2063	endif
2064	endif
2065
2066	if(iteM>=MAX_FIND_ITER) then
2067	! C_%stable_da=.FALSE.
2068	! IF(iteM==MAX_FIND_ITER+100) THEN
2069	! write(6,*) " Unstable in find_orbit without TPSA"
2070	messagelost= "Maximum number of iterations in find_orbit without TPSA"
2071	xlost=fix
2072	check_stable=my_false
2073	! ENDIF
2074	ITE=0
2075	endif
2076
2077	if(ite.eq.1) then
2078	GOTO 3
2079
2080	endif
2081	c_%APERTURE_FLAG=APERTURE
2082
2083	END SUBROUTINE FIND_ORBIT_LAYOUT_noda
2084
2085
2086	SUBROUTINE fit_all_bends(r,state)
2087	IMPLICIT NONE
2088	TYPE(layout),target,INTENT(INOUT):: r
2089	TYPE(internal_state), intent(in):: state
2090	type(fibre),pointer :: p
2091	integer i
2092
2093	p=>r%start
2094
2095	do i=1,r%n
2096	if(p%mag%p%b0/=0.0_dp) call fit_bare_bend(p,state)
2097	p=>p%next
2098	enddo
2099
2100	end SUBROUTINE fit_all_bends
2101
2102	SUBROUTINE fit_bare_bend(f,state,next)
2103	IMPLICIT NONE
2104	TYPE(fibre),INTENT(INOUT):: f
2105	TYPE(real_8) y(6)
2106	TYPE(internal_state), intent(in):: state
2107	! integer,optional,target :: charge
2108	real(dp) kf,x(6),xdix,xdix0,tiny
2109	integer ite
2110	logical(lp), optional :: next
2111	logical(lp) nex
2112	nex=my_false
2113	if(present(next)) nex=next
2114	tiny=1e-40_dp
2115	xdix0=1e4_dp*DEPS_tracking
2116
2117	KF=0.0_dp ;
2118	F%MAGP%BN(1)%KIND=3
2119	F%MAGP%BN(1)%I=1
2120	if(nex) then
2121	F%next%MAGP%BN(1)%KIND=3
2122	F%next%MAGP%BN(1)%I=1
2123	endif
2124
2125	CALL INIT(1,1)
2126
2127	CALL ALLOC(Y)
2128
2129	3 continue
2130	X=0.0_dp
2131	Y=X
2132	CALL TRACK(f,Y,+state) !,CHARGE)
2133	if(nex) CALL TRACK(f%next,Y,+state) !,CHARGE)
2134	x=y
2135	! write(6,'(A10,6(1X,G14.7))') " ORBIT IS ",x
2136	kf=-(y(2).sub.'0')/(y(2).sub.'1')
2137	xdix=abs(y(2).sub.'0')
2138	f%MAGP%BN(1) = f%MAGP%BN(1)+KF
2139	f%MAG%BN(1) = f%MAG%BN(1)+KF
2140	if(nex) then
2141	f%next%MAGP%BN(1) = f%next%MAGP%BN(1)+KF
2142	f%next%MAG%BN(1) = f%next%MAG%BN(1)+KF
2143	endif
2144
2145	CALL ADD(f,1,1,0.0_dp) !etienne
2146	if(nex) CALL ADD(f%next,1,1,0.0_dp) !etienne
2147
2148	if(xdix.gt.deps_tracking) then
2149	ite=1
2150	else
2151	if(xdix.ge.xdix0.or.xdix<=tiny) then
2152	ite=0
2153	else
2154	ite=1
2155	xdix0=xdix
2156	endif
2157	endif
2158	if(ite.eq.1) GOTO 3
2159
2160	F%MAGP%BN(1)%KIND=1
2161	F%MAGP%BN(1)%I=0
2162	if(nex) then
2163	F%next%MAGP%BN(1)%KIND=1
2164	F%next%MAGP%BN(1)%I=0
2165	endif
2166	! write(6,'(A10,1(1X,g14.7))') " BN(1) IS ", f%MAG%BN(1)
2167
2168
2169	CALL KILL(Y)
2170
2171
2172	end SUBROUTINE fit_bare_bend
2173
2174	SUBROUTINE track_aperture(r,my_state,beta,dbeta,tuneold,ib,ITMAX,emit0,aper,pos,nturn,FILENAME,FILEtune,FILESMEAR,resmax)
2175	IMPLICIT NONE
2176	TYPE(INTERNAL_STATE), intent(IN):: my_STATE
2177	TYPE(INTERNAL_STATE) STATE
2178	TYPE(layout),target, intent(inout) :: R
2179	REAL(DP), pointer :: BETA(:,:,:)
2180	integer pos,nturn,i,flag,ib,MF,mft,j,resmax,it,I1,no
2181	real(dp) closed(6),MAT(6,6),AI(6,6),A(6,6),emit(2),emit0(6),aper(2),x(6),xn(6),dbeta,tuneold(:)
2182	real(dp) ra(2),tunenew(2),xda(lnv)
2183	CHARACTER(*) FILENAME,FILEtune,FILESMEAR
2184	real(dp), allocatable :: dat(:,:),dats(:,:),SMEAR(:,:)
2185	REAL(DP) JMin(2),JMAX(2), tune1(2),tune2(2),tot_tune(2),epsi,scas(2),scau,scat(2)
2186	integer itmax
2187	type(damap) id
2188	type(tree) monkey
2189
2190
2191
2192	epsi=1.0_dp/nturn
2193	STATE=((((my_state+nocavity0)+delta0)+only_4d0)-RADIATION0)
2194	allocate(dat(0:nturn,6),dats(0:nturn,6))
2195	allocate(SMEAR(ITMAX,8))
2196	CLOSED=0.0_dp
2197
2198	call FILL_BETA(r,my_state,pos,BETA,IB,DBETA,tuneold,tunenew,a,ai,mat,closed)
2199	write(6,) " ****************************************************************"
2200	write(6,*) " Tracking with Normalized Aperture "
2201	write(6,*) " Tunes = ",tunenew(1:2)
2202	if(pos==0) then
2203	write(6,*) " give no "
2204	read(5,*) no
2205	call init(no,2,0,0)
2206	call alloc(id); call alloc(monkey)
2207	call kanalnummer(mf)
2208	open(unit=mf,file='map.dat')
2209	call read(id,mf)
2210	id=closed
2211	monkey=id
2212	call kill(id)
2213	close(mf)
2214	xda=0.0_dp
2215	endif
2216	scau=1.0_dp
2217	scas=0.0_dp
2218	dats=0.0_dp
2219	SMEAR=0.0_dp
2220	it=0
2221	CALL KANALNUMMER(MFt)
2222	OPEN(UNIT=MFt,FILE=FILEtune)
2223	1001 continue
2224	it=it+1
2225	write(6,*) " iteration ",it
2226	IF(IT==ITMAX+1) GOTO 1002
2227
2228	scat(1)=(emit0(1)+ it*emit0(3))/aper(1)
2229	scat(2)=(emit0(2)+ it*emit0(6))/aper(2)
2230	! scat=(scau+scas)/two ! etienne
2231	dat=0.0_dp
2232
2233	xn=0.0_dp
2234	JMAX=0.0_dp
2235	JMIN=mybig
2236	emit=scat*aper
2237	write(6,*) " %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%"
2238	write(6,*) " Initial emit = ", emit(1:2)," scale = ",scat
2239	write(6,*) " %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%"
2240	xn(2)=sqrt(emit(1))
2241	xn(4)=sqrt(emit(2))
2242	X=0.0_dp
2243
2244	x=matmul(a,xn)+closed
2245
2246	WRITE(6,*) " INITIAL X"
2247	WRITE(6,200) X
2248	! x=zero
2249	! x(1)=xn(2)/sqrt(a(2,1)2+a(2,2)2)
2250	! x(3)=xn(4)/sqrt(a(4,3)2+a(4,4)2)
2251	! Etienne
2252
2253	dat(0,1:4)=xn(1:4)
2254	ra(1)=xn(1)2+xn(2)2
2255	ra(2)=xn(3)2+xn(4)2
2256	dat(0,5:6)=ra(1:2)
2257
2258	flag=0
2259	do i=1,nturn
2260
2261	if(pos/=0) then
2262	flag=track_flag(r,x,pos,state)
2263	else
2264	flag=0
2265	xda(1:4)=x(1:4)
2266	xda=monkey*xda
2267	x(1:4)=xda(1:4)
2268	endif
2269	write(80,*) x(1:4)
2270	if(flag/=0) exit
2271	xn=matmul(ai,(x-closed))
2272	ra(1)=xn(1)2+xn(2)2
2273	ra(2)=xn(3)2+xn(4)2
2274	IF(RA(1)>JMAX(1)) JMAX(1)=RA(1)
2275	IF(RA(2)>JMAX(2)) JMAX(2)=RA(2)
2276	IF(RA(1)<JMIN(1)) JMIN(1)=RA(1)
2277	IF(RA(2)<JMIN(2)) JMIN(2)=RA(2)
2278	if(ra(1)>aper(1)) then
2279	flag=101
2280	exit
2281	endif
2282	if(ra(2)>aper(2)) then
2283	flag=102
2284	exit
2285	endif
2286	dat(i,1:4)=xn(1:4)
2287	dat(i,5:6)=ra(1:2)
2288	enddo
2289	WRITE(6,*) " MAXIMUM RADIUS SQUARE = ",JMAX
2290	WRITE(6,*) " MAXIMUM/INITIAL = ",JMAX(1:2)/emit(1:2)
2291	WRITE(6,*) " MINIMUM RADIUS SQUARE = ",JMIN
2292	WRITE(6,*) " MINIMUM/INITIAL = ",JMIN(1:2)/emit(1:2)
2293	WRITE(6,) " SMEAR = ",2.0_dp(JMAX-JMIN)/(JMAX+JMIN)
2294	SMEAR(IT,1:2)=EMIT(1:2)
2295	SMEAR(IT,3:4)=JMIN(1:2)
2296	SMEAR(IT,5:6)=JMAX(1:2)
2297	SMEAR(IT,7:8)=2.0_dp*(JMAX-JMIN)/(JMAX+JMIN)
2298	if(flag/=0)then
2299	! scau=scat
2300	IF(fLAG==101) THEN
2301	write(6,*) " UNSTABLE AT X-NORMALIZED APERTURE "
2302	! write(mf,*) "UNSTABLE AT X-NORMALIZED APERTURE "
2303	ELSEIF(FLAG==102) THEN
2304	! write(mf,*) "UNSTABLE AT Y-NORMALIZED APERTURE "
2305	write(6,*) " UNSTABLE AT Y-NORMALIZED APERTURE "
2306	ELSE
2307	! write(mf,*) "UNSTABLE: DYNAMIC APERTURE "
2308	write(6,*) " UNSTABLE: DYNAMIC APERTURE "
2309	ENDIF
2310	goto 1002 ! etienne
2311	if(abs(scau-scas(1))<=emit0(3)) then
2312	goto 1002
2313	else
2314	goto 1001
2315	endif
2316	ELSE
2317	write(6,*) " STABLE "
2318
2319	! write(6,*) "tunes of the ray " ,tot_tune(1:2)
2320	WRITE(6,201) EMIT,APER, TUNEnew(1:2),DBETA
2321
2322	scas=scat
2323	dats=dat
2324
2325	! RESONANCE
2326	tot_tune=0.0_dp
2327	xn(1:4)=dats(0,1:4)
2328	tune1(1)=atan2(-xn(2),xn(1))/twopi
2329	tune1(2)=atan2(-xn(4),xn(3))/twopi
2330	if(tune1(1)<0.0_dp) tune1(1)=tune1(1)+1.0_dp
2331	if(tune1(2)<0.0_dp) tune1(2)=tune1(2)+1.0_dp
2332	DO I1=0,NTURN
2333	xn(1:4)=dats(i1,1:4)
2334	tune2(1)=atan2(-xn(2),xn(1))/twopi
2335	tune2(2)=atan2(-xn(4),xn(3))/twopi
2336	if(tune2(1)<0.0_dp) tune2(1)=tune2(1)+1.0_dp
2337	if(tune2(2)<0.0_dp) tune2(2)=tune2(2)+1.0_dp
2338	tune1=tune2-tune1
2339	if(tune1(1)<0.0_dp) tune1(1)=tune1(1)+1.0_dp
2340	if(tune1(2)<0.0_dp) tune1(2)=tune1(2)+1.0_dp
2341	tot_tune =tot_tune+tune1
2342	tune1=tune2
2343	ENDDO
2344	tot_tune=tot_tune/nturn
2345	do i1=0,resmax
2346	do j=-resmax,resmax
2347	dbeta=i1tot_tune(1)+jtot_tune(2)
2348	if(abs(dbeta-nint(dbeta))<epsi ) then
2349	if(i1+j/=0) write(6,*) i1,j,dbeta," <--- here "
2350	else
2351	if(i1+j/=0) write(6,*)i1,j,dbeta
2352	endif
2353	enddo
2354	enddo
2355	write(mft,*) " emit = ",emit(1:2)
2356	write(mft,*) " tunes = ",tot_tune(1:2)
2357	do i=0,resmax
2358	do j=-resmax,resmax
2359	dbeta=itot_tune(1)+jtot_tune(2)
2360	if(abs(dbeta-nint(dbeta))<epsi ) then
2361	if(i+j/=0) write(mft,*) i,j,dbeta," <--- here "
2362	else
2363	if(i+j/=0) write(mft,*)i,j,dbeta
2364	endif
2365	enddo
2366	enddo
2367	! PAUSE 100
2368	! END RESONANCE
2369
2370
2371
2372	if(abs(scau-scas(1))<=emit0(3)) then
2373	goto 1002
2374	else
2375	goto 1001
2376	endif
2377
2378	ENDIF
2379
2380
2381	1002 continue
2382	CALL KANALNUMMER(MF)
2383
2384	OPEN(UNIT=MF,FILE=FILENAME)
2385
2386	WRITE(MF,201) EMIT,APER, TUNEnew(1:2),DBETA
2387	tot_tune=0.0_dp
2388	xn(1:4)=dats(0,1:4)
2389	tune1(1)=atan2(-xn(2),xn(1))/twopi
2390	tune1(2)=atan2(-xn(4),xn(3))/twopi
2391	if(tune1(1)<0.0_dp) tune1(1)=tune1(1)+1.0_dp
2392	if(tune1(2)<0.0_dp) tune1(2)=tune1(2)+1.0_dp
2393
2394	DO I=0,NTURN
2395	WRITE(MF,200)DATs(I,1:6)
2396	xn(1:4)=dats(i,1:4)
2397	tune2(1)=atan2(-xn(2),xn(1))/twopi
2398	tune2(2)=atan2(-xn(4),xn(3))/twopi
2399	if(tune2(1)<0.0_dp) tune2(1)=tune2(1)+1.0_dp
2400	if(tune2(2)<0.0_dp) tune2(2)=tune2(2)+1.0_dp
2401	tune1=tune2-tune1
2402	if(tune1(1)<0.0_dp) tune1(1)=tune1(1)+1.0_dp
2403	if(tune1(2)<0.0_dp) tune1(2)=tune1(2)+1.0_dp
2404	tot_tune =tot_tune+tune1
2405	tune1=tune2
2406	ENDDO
2407	tot_tune=tot_tune/nturn
2408	CLOSE(MF)
2409
2410	CLOSE(mft)
2411
2412	CALL KANALNUMMER(MF)
2413	OPEN(UNIT=MF,FILE=FILESMEAR)
2414	WRITE(MF,) " ITERATION EMIT0(1:2) JMIN(1:2) JMAX(1:2) 2.0_dp(JMAX-JMIN)/(JMAX+JMIN)"
2415	DO I=1,ITMAX
2416	WRITE(MF,202)I, SMEAR(I,1:8)
2417	ENDDO
2418
2419	CLOSE(mf)
2420
2421	emit0(1:2) =scas*aper
2422	emit0(4:5)=tunenew(1:2)
2423
2424	202 FORMAT(1X,I4,8(1X,D18.11))
2425	201 FORMAT(9(1X,D18.11))
2426	200 FORMAT(6(1X,D18.11))
2427	deallocate(dat,dats,SMEAR)
2428	WRITE(6,*) " MAXIMUM RADIUS SQUARE = ",JMAX
2429	WRITE(6,*) " MAXIMUM/INITIAL = ",JMAX(1:2)/emit(1:2)
2430	WRITE(6,*) " MINIMUM RADIUS SQUARE = ",JMIN
2431	WRITE(6,*) " MINIMUM/INITIAL = ",JMIN(1:2)/emit(1:2)
2432	WRITE(6,) " SMEAR = ",2.0_dp(JMAX-JMIN)/(JMAX+JMIN)
2433	write(6,) " ****************************************************************"
2434	if(pos==0) call kill(monkey)
2435	end SUBROUTINE track_aperture
2436
2437
2438	SUBROUTINE THIN_LENS_resplit(R,THIN,even,lim,lmax0,xbend,sexr,fib,useknob) ! A re-splitting routine
2439	IMPLICIT NONE
2440	INTEGER NTE
2441	TYPE(layout),target, intent(inout) :: R
2442	real(dp), OPTIONAL, intent(inout) :: THIN
2443	real(dp), OPTIONAL, intent(in) :: lmax0
2444	real(dp), OPTIONAL, intent(in) ::xbend
2445	real(dp), OPTIONAL, intent(in) ::sexr
2446	type(fibre), OPTIONAL, target :: fib
2447	logical(lp), OPTIONAL :: useknob
2448	real(dp) gg,RHOI,XL,QUAD,THI,lm,dl,ggbt,xbend1,gf(7),sexr0,quad0
2449	INTEGER M1,M2,M3, MK1,MK2,MK3,limit(2),parity,inc,nst_tot,ntec,ii,metb,sexk
2450	integer incold ,parityold
2451	integer, optional :: lim(2)
2452	logical(lp) MANUAL,eject,doit,DOBEND
2453	TYPE (fibre), POINTER :: C
2454	logical(lp),optional :: even
2455	type(layout), pointer :: L
2456	logical f1,f2
2457
2458	sexr0=0.0_dp
2459
2460
2461	if(present(sexr)) sexr0=sexr
2462
2463	f1=.false.
2464	f2=.false.
2465
2466	if(associated(r%parent_universe)) then
2467	f1=.true.
2468	l=>r%parent_universe%start
2469	do ii=1,r%parent_universe%n
2470	call kill(l%t)
2471	l=>l%next
2472	enddo
2473	elseif(associated(r%t)) then
2474	call kill(r%t)
2475	f2=.true.
2476	endif
2477	! logical(lp) doneit
2478	nullify(C)
2479	parity=0
2480	inc=0
2481	lm=1.0e38_dp
2482	ntec=0
2483	max_ds=0.0_dp
2484	xbend1=-1.0_dp
2485
2486	if(present(xbend)) xbend1=xbend
2487	if(present(lmax0)) lm=abs(lmax0)
2488	if(present(even)) then
2489	inc=1
2490	if(even) then
2491	parity=0
2492	else
2493	parity=1
2494	endIf
2495	endif
2496	parityold=parity
2497	incold=inc
2498	! CALL LINE_L(R,doneit)
2499
2500	MANUAL=.FALSE.
2501	eject=.FALSE.
2502
2503	THI=R%THIN
2504	IF(PRESENT(THIN)) THI=THIN
2505
2506	IF(THI<=0) MANUAL=.TRUE.
2507
2508
2509	IF(MANUAL) THEN
2510	write(6,*) "thi: thin lens factor (THI<0 TO STOP), sextupole factor and Bend factor "
2511	read(5,*) thi,sexr0,xbend1
2512	IF(THI<0) eject=.true.
2513	ENDIF
2514
2515	1001 CONTINUE
2516
2517	limit(1)=4
2518	limit(2)=18
2519	if(present(lim)) limit=lim
2520	if(sixtrack_compatible) then
2521	limit(1)=1000000
2522	limit(2)=1000001
2523	endif
2524	! limit0(1)=limit(1)
2525	! limit0(2)=limit(2)
2526
2527	M1=0
2528	M2=0
2529	M3=0
2530	MK1=0
2531	MK2=0
2532	MK3=0
2533	r%NTHIN=0
2534	nst_tot=0
2535	sexk=0
2536	C=>R%START
2537	do ii=1,r%n ! WHILE(ASSOCIATED(C))
2538	doit=(C%MAG%KIND==kind1.or.C%MAG%KIND==kind2.or.C%MAG%KIND==kind4.or.C%MAG%KIND==kind5)
2539	doit=DOIT.OR.(C%MAG%KIND==kind6.or.C%MAG%KIND==kind7)
2540	DOIT=DOIT.OR.(C%MAG%KIND==kind10.or.C%MAG%KIND==kind16)
2541	DOIT=DOIT.OR.(C%MAG%KIND==kind17)
2542	doit=doit.and.C%MAG%recut
2543
2544	if(doit) then
2545	select case(C%MAG%P%METHOD)
2546	CASE(2)
2547	M1=M1+1
2548	MK1=MK1+C%MAG%P%NST
2549	CASE(4)
2550	M2=M2+1
2551	MK2=MK2+3*C%MAG%P%NST
2552	CASE(6)
2553	M3=M3+1
2554	MK3=MK3+7*C%MAG%P%NST
2555	END SELECT
2556	r%NTHIN=r%NTHIN+1 !C%MAG%NST
2557	endif
2558	NST_tot=NST_tot+C%MAG%P%nst
2559
2560	C=>C%NEXT
2561
2562	enddo
2563	write(6,*) "Previous of cutable Elements ",r%NTHIN
2564	write(6,*) "METHOD 2 ",M1,MK1
2565	write(6,*) "METHOD 4 ",M2,MK2
2566	write(6,*) "METHOD 6 ",M3,MK3
2567	write(6,*) "number of Slices ", MK1+MK2+MK3
2568	write(6,*) "Total NST ", NST_tot
2569
2570	if(eject) then
2571	! limit(1)=limit0(1)
2572	! limit(2)=limit0(2)
2573	return
2574	endif
2575	M1=0
2576	M2=0
2577	M3=0
2578	MK1=0
2579	MK2=0
2580	MK3=0
2581	! CAVITY FOCUSING
2582	! TEAPOT SPLITTING....
2583	ggbt=0.0_dp
2584	r%NTHIN=0
2585	r%THIN=THI
2586
2587	nst_tot=0
2588	C=>R%START
2589	do ii=1,r%n ! WHILE(ASSOCIATED(C))
2590	doit=.true.
2591	if(present(useknob)) then
2592	if(useknob) then
2593	doit=c%magp%knob
2594	else
2595	doit=.not.c%magp%knob
2596	endif
2597	endif
2598	if(present(fib)) then
2599	doit=associated(c,fib)
2600	endif
2601	if(.not.doit) then
2602	NST_tot=NST_tot+C%MAG%P%nst
2603	if(c%mag%even) then
2604	parity=parityold
2605	inc=incold
2606	endif
2607	C=>C%NEXT
2608	cycle
2609	endif
2610
2611	if(c%mag%even) then
2612	parity=0
2613	inc=1
2614	endif
2615
2616
2617	! if(doit) then
2618
2619	select case(resplit_cutting)
2620
2621	case(0)
2622
2623	doit=(C%MAG%KIND==kind1.or.C%MAG%KIND==kind2.or.C%MAG%KIND==kind4.or.C%MAG%KIND==kind5)
2624	doit=DOIT.OR.(C%MAG%KIND==kind6.or.C%MAG%KIND==kind7)
2625	DOIT=DOIT.OR.(C%MAG%KIND==kind10.or.C%MAG%KIND==kind16)
2626	DOIT=DOIT.OR.(C%MAG%KIND==kind17)
2627	doit=doit.and.C%MAG%recut
2628	if(doit) then
2629	xl=C%MAG%L
2630	RHOI=0.0_dp
2631	QUAD=0.0_dp
2632	QUAD0=0.0_dp
2633	IF(C%MAG%P%NMUL>=1) THEN
2634	! RHOI=C%MAG%P%B0
2635	RHOI=abs(C%MAG%bn(1))+abs(C%MAG%an(1))
2636	endif
2637	IF(C%MAG%P%NMUL>=2) THEN
2638	QUAD=SQRT(C%MAG%BN(2)2+C%MAG%AN(2)2)
2639	IF(C%MAG%P%NMUL>=3) THEN
2640	quad0=SQRT(C%MAG%BN(3)2+C%MAG%AN(3)2)*sexr0
2641	QUAD=QUAD+quad0
2642	endif
2643	ELSE
2644	QUAD=0.0_dp
2645	ENDIF
2646	if(C%MAG%KIND==kind5.or.C%MAG%KIND==kind17) then
2647	quad=quad+(C%MAG%b_sol)**2/4.0_dp+abs(C%MAG%b_sol/2.0_dp)
2648	endif
2649
2650	DOBEND=MY_FALSE
2651	IF(xbend1>0.0_dp) THEN
2652	IF(C%MAG%KIND==kind10) THEN
2653	IF(C%MAG%TP10%DRIFTKICK) THEN
2654	DOBEND=MY_TRUE
2655	ENDIF
2656	ENDIF
2657	IF(C%MAG%KIND==kind16.OR.C%MAG%KIND==kind20) THEN
2658	IF(C%MAG%K16%DRIFTKICK) THEN
2659	DOBEND=MY_TRUE
2660	ENDIF
2661	ENDIF
2662	if(rhoi/=0.0_dp.and.radiation_bend_split)DOBEND=MY_TRUE
2663	ENDIF
2664	! ETIENNE
2665	GG=XL(RHOI*2+ABS(QUAD))
2666	GG=GG/THI
2667	NTE=INT(GG)
2668	sexk=sexk+xl*ABS(QUAD0)/thi
2669	metb=0
2670	if(dobend) then
2671	call check_bend(xl,gg,rhoi,xbend1,gf,metb)
2672	if(gf(metb)>gg) then
2673	gg=gf(metb)
2674	NTE=INT(GG)
2675	ggbt=ggbt+NTE
2676	else
2677	metb=0
2678	endif
2679	endif
2680
2681
2682	IF(NTE.LT.limit(1).or.metb==2) THEN
2683	M1=M1+1
2684	IF(NTE.EQ.0) NTE=1
2685	if(mod(nte,2)/=parity) nte=nte+inc
2686	C%MAG%P%NST=NTE
2687	C%MAG%P%METHOD=2
2688	MK1=MK1+NTE
2689	ELSEIF((NTE.GE.limit(1).AND.NTE.LT.limit(2)).or.metb==4) THEN
2690	M2=M2+1
2691	NTE=NTE/3
2692	IF(NTE.EQ.0) NTE=1
2693	if(mod(nte,2)/=parity) nte=nte+inc
2694	C%MAG%P%NST=NTE
2695	C%MAG%P%METHOD=4
2696	MK2=MK2+NTE*3
2697	ELSEIF(NTE.GE.limit(2).or.metb==6) THEN
2698	M3=M3+1
2699	NTE=NTE/7
2700	IF(NTE.EQ.0) NTE=1
2701	if(mod(nte,2)/=parity) nte=nte+inc
2702	C%MAG%P%NST=NTE
2703	C%MAG%P%METHOD=6
2704	MK3=MK3+NTE*7
2705	ENDIF
2706
2707	r%NTHIN=r%NTHIN+1 !C%MAG%NST
2708	! write(6,*)"nte>ntec", nte,ntec
2709	call add(C%MAG,C%MAG%P%nmul,1,0.0_dp)
2710	call COPY(C%MAG,C%MAGP)
2711	if(gg>0.0_dp) then
2712	if(c%mag%l/c%mag%p%nst>max_ds) max_ds=c%mag%l/c%mag%p%nst
2713	endif
2714
2715	endif ! doit
2716
2717	case(1)
2718
2719	doit=(C%MAG%KIND==kind1.or.C%MAG%KIND==kind2.or.C%MAG%KIND==kind4.or.C%MAG%KIND==kind5)
2720	doit=DOIT.OR.(C%MAG%KIND==kind6.or.C%MAG%KIND==kind7)
2721	DOIT=DOIT.OR.(C%MAG%KIND==kind10.or.C%MAG%KIND==kind16)
2722	DOIT=DOIT.OR.(C%MAG%KIND==kind17)
2723	doit=doit.and.C%MAG%recut
2724
2725	if(doit) then
2726	xl=C%MAG%L
2727	RHOI=0.0_dp
2728	QUAD=0.0_dp
2729	QUAD0=0.0_dp
2730	IF(C%MAG%P%NMUL>=1) THEN
2731	! RHOI=C%MAG%P%B0
2732	RHOI=abs(C%MAG%bn(1))+abs(C%MAG%an(1))
2733	endif
2734	IF(C%MAG%P%NMUL>=2) THEN
2735	QUAD=SQRT(C%MAG%BN(2)2+C%MAG%AN(2)2)
2736	IF(C%MAG%P%NMUL>=3) THEN
2737	quad0=SQRT(C%MAG%BN(3)2+C%MAG%AN(3)2)*sexr0
2738	QUAD=QUAD+quad0
2739	endif
2740	ELSE
2741	QUAD=0.0_dp
2742	ENDIF
2743	if(C%MAG%KIND==kind5.or.C%MAG%KIND==kind17) then
2744	quad=quad+(C%MAG%b_sol)**2/4.0_dp+abs(C%MAG%b_sol/2.0_dp)
2745	endif
2746
2747	DOBEND=MY_FALSE
2748	IF(xbend1>0.0_dp) THEN
2749	IF(C%MAG%KIND==kind10) THEN
2750	IF(C%MAG%TP10%DRIFTKICK) THEN
2751	DOBEND=MY_TRUE
2752	ENDIF
2753	ENDIF
2754	IF(C%MAG%KIND==kind16.OR.C%MAG%KIND==kind20) THEN
2755	IF(C%MAG%K16%DRIFTKICK) THEN
2756	DOBEND=MY_TRUE
2757	ENDIF
2758	ENDIF
2759	if(rhoi/=0.0_dp.and.radiation_bend_split)DOBEND=MY_TRUE
2760	ENDIF
2761	! ETIENNE
2762	GG=XL(RHOI*2+ABS(QUAD))
2763	GG=GG/THI
2764	sexk=sexk+xl*ABS(QUAD0)/thi
2765	NTE=INT(GG)
2766	metb=0
2767	if(dobend) then
2768	call check_bend(xl,gg,rhoi,xbend1,gf,metb)
2769	if(gf(metb)>gg) then
2770	gg=gf(metb)
2771	NTE=INT(GG)
2772	ggbt=ggbt+NTE
2773	else
2774	metb=0
2775	endif
2776	endif
2777
2778
2779	IF(NTE.LT.limit(1).or.metb==2) THEN
2780	M1=M1+1
2781	IF(NTE.EQ.0) NTE=1
2782	if(mod(nte,2)/=parity) nte=nte+inc
2783	C%MAG%P%NST=NTE
2784	C%MAG%P%METHOD=2
2785	MK1=MK1+NTE
2786	ELSEIF((NTE.GE.limit(1).AND.NTE.LT.limit(2)).or.metb==4) THEN
2787	M2=M2+1
2788	NTE=NTE/3
2789	IF(NTE.EQ.0) NTE=1
2790	if(mod(nte,2)/=parity) nte=nte+inc
2791	C%MAG%P%NST=NTE
2792	C%MAG%P%METHOD=4
2793	MK2=MK2+NTE*3
2794	ELSEIF(NTE.GE.limit(2).or.metb==6) THEN
2795	M3=M3+1
2796	NTE=NTE/7
2797	IF(NTE.EQ.0) NTE=1
2798	if(mod(nte,2)/=parity) nte=nte+inc
2799	C%MAG%P%NST=NTE
2800	C%MAG%P%METHOD=6
2801	MK3=MK3+NTE*7
2802	ENDIF
2803
2804
2805	r%NTHIN=r%NTHIN+1 !C%MAG%NST
2806
2807	if(present(lmax0).and.c%mag%kind==kind1) then
2808	dl=(C%MAG%P%ld/C%MAG%P%nst)
2809	if(dl>lm*fuzzy_split) then
2810	ntec=int(C%MAG%P%ld/lm)+1
2811	if(mod(ntec,2)/=parity) ntec=ntec+inc
2812	C%MAG%P%NST=ntec
2813	endif
2814	endif
2815
2816	call add(C%MAG,C%MAG%P%nmul,1,0.0_dp)
2817	call COPY(C%MAG,C%MAGP)
2818	! if(gg>zero) then
2819	if(c%mag%l/c%mag%p%nst>max_ds) max_ds=c%mag%l/c%mag%p%nst
2820	! endif
2821	endif
2822
2823	case(2)
2824
2825
2826	doit=(C%MAG%KIND==kind1.or.C%MAG%KIND==kind2.or.C%MAG%KIND==kind4.or.C%MAG%KIND==kind5)
2827	doit=DOIT.OR.(C%MAG%KIND==kind6.or.C%MAG%KIND==kind7)
2828	DOIT=DOIT.OR.(C%MAG%KIND==kind10.or.C%MAG%KIND==kind16)
2829	DOIT=DOIT.OR.(C%MAG%KIND==kind17)
2830	doit=doit.and.C%MAG%recut
2831
2832	if(doit) then
2833	xl=C%MAG%L
2834	RHOI=0.0_dp
2835	QUAD=0.0_dp
2836	QUAD0=0.0_dp
2837	IF(C%MAG%P%NMUL>=1) THEN
2838	! RHOI=C%MAG%P%B0
2839	RHOI=abs(C%MAG%bn(1))+abs(C%MAG%an(1))
2840	endif
2841	IF(C%MAG%P%NMUL>=2) THEN
2842	QUAD=SQRT(C%MAG%BN(2)2+C%MAG%AN(2)2)
2843	IF(C%MAG%P%NMUL>=3) THEN
2844	quad0=SQRT(C%MAG%BN(3)2+C%MAG%AN(3)2)*sexr0
2845	QUAD=QUAD+quad0
2846	endif
2847	ELSE
2848	QUAD=0.0_dp
2849	ENDIF
2850	if(C%MAG%KIND==kind5.or.C%MAG%KIND==kind17) then
2851	quad=quad+(C%MAG%b_sol)**2/4.0_dp+abs(C%MAG%b_sol/2.0_dp)
2852	endif
2853
2854	DOBEND=MY_FALSE
2855	IF(xbend1>0.0_dp) THEN
2856	IF(C%MAG%KIND==kind10) THEN
2857	IF(C%MAG%TP10%DRIFTKICK) THEN
2858	DOBEND=MY_TRUE
2859	ENDIF
2860	ENDIF
2861	IF(C%MAG%KIND==kind16.OR.C%MAG%KIND==kind20) THEN
2862	IF(C%MAG%K16%DRIFTKICK) THEN
2863	DOBEND=MY_TRUE
2864	ENDIF
2865	ENDIF
2866	if(rhoi/=0.0_dp.and.radiation_bend_split)DOBEND=MY_TRUE
2867	ENDIF
2868	! ETIENNE
2869	GG=XL(RHOI*2+ABS(QUAD))
2870	GG=GG/THI
2871	NTE=INT(GG)
2872	sexk=sexk+xl*ABS(QUAD0)/thi
2873	metb=0
2874	if(dobend) then
2875	call check_bend(xl,gg,rhoi,xbend1,gf,metb)
2876	if(gf(metb)>gg) then
2877	gg=gf(metb)
2878	NTE=INT(GG)
2879	ggbt=ggbt+NTE
2880	else
2881	metb=0
2882	endif
2883	endif
2884
2885
2886	IF(NTE.LT.limit(1).or.metb==2) THEN
2887	M1=M1+1
2888	IF(NTE.EQ.0) NTE=1
2889	if(mod(nte,2)/=parity) nte=nte+inc
2890	C%MAG%P%NST=NTE
2891	C%MAG%P%METHOD=2
2892	MK1=MK1+NTE
2893	ELSEIF((NTE.GE.limit(1).AND.NTE.LT.limit(2)).or.metb==4) THEN
2894	M2=M2+1
2895	NTE=NTE/3
2896	IF(NTE.EQ.0) NTE=1
2897	if(mod(nte,2)/=parity) nte=nte+inc
2898	C%MAG%P%NST=NTE
2899	C%MAG%P%METHOD=4
2900	MK2=MK2+NTE*3
2901	ELSEIF(NTE.GE.limit(2).or.metb==6) THEN
2902	M3=M3+1
2903	NTE=NTE/7
2904	IF(NTE.EQ.0) NTE=1
2905	if(mod(nte,2)/=parity) nte=nte+inc
2906	C%MAG%P%NST=NTE
2907	C%MAG%P%METHOD=6
2908	MK3=MK3+NTE*7
2909	ENDIF
2910
2911
2912	r%NTHIN=r%NTHIN+1 !C%MAG%NST
2913	! write(6,*)"nte>ntec", nte,ntec
2914	if(nte>ntec.or.(.not.present(lmax0)) ) then
2915	call add(C%MAG,C%MAG%P%nmul,1,0.0_dp)
2916	call COPY(C%MAG,C%MAGP)
2917	endif
2918	! if(gg>zero) then
2919	! if(c%mag%l/c%mag%p%nst>max_ds) max_ds=c%mag%l/c%mag%p%nst
2920	! endif
2921
2922	if(present(lmax0)) then
2923	dl=(C%MAG%P%ld/C%MAG%P%nst)
2924	if(dl>lm*fuzzy_split.and.C%MAG%KIND/=kindpa) then
2925	nte=int(C%MAG%P%ld/lm)+1
2926	if(mod(nte,2)/=parity) nte=nte+inc
2927	if(nte > C%MAG%P%NST ) then
2928	C%MAG%P%NST=nte
2929	call add(C%MAG,C%MAG%P%nmul,1,0.0_dp)
2930	call COPY(C%MAG,C%MAGP)
2931	endif
2932
2933	elseif(dl>lm.and.C%MAG%KIND==kindpa) then
2934	write(6,*) " Pancake cannot be recut "
2935	endif
2936	endif
2937	if(c%mag%l/c%mag%p%nst>max_ds) max_ds=c%mag%l/c%mag%p%nst
2938
2939
2940	endif
2941
2942	case default
2943	stop 988
2944	end select
2945
2946
2947	! endif
2948	NST_tot=NST_tot+C%MAG%P%nst
2949	if(c%mag%even) then
2950	parity=parityold
2951	inc=incold
2952	endif
2953	C=>C%NEXT
2954
2955	enddo ! end of do WHILE
2956
2957
2958	write(6,*) "Present of cutable Elements ",r%NTHIN
2959	write(6,*) "METHOD 2 ",M1,MK1
2960	write(6,*) "METHOD 4 ",M2,MK2
2961	write(6,*) "METHOD 6 ",M3,MK3
2962	write(6,*) "number of Slices ", MK1+MK2+MK3
2963	write(6,*) "Total NST ", NST_tot
2964	if(radiation_bend_split) then
2965	write(6,*) "Total NST due to Bend Closed Orbit ", int(ggbt)
2966	write(6,*) "Restricted to method=2 for radiation or spin "
2967	else
2968	write(6,*) "Total NST due to Bend Closed Orbit ", int(ggbt)
2969	endif
2970	write(6,*) "Total NST due to Sextupoles ", sexk
2971	write(6,*) "Biggest ds ", max_ds
2972
2973
2974
2975	IF(MANUAL) THEN
2976	write(6,*) "thi: thin lens factor (THI<0 TO STOP), sextupole factor and Bend factor "
2977	read(5,*) thi,sexr0, xbend1
2978	IF(THI<0) THEN
2979	THI=R%THIN
2980	! limit(1)=limit0(1)
2981	! limit(2)=limit0(2)
2982	if(f1) then
2983	l=>r%parent_universe%start
2984	do ii=1,r%parent_universe%n
2985	call make_node_layout(l)
2986	l=>l%next
2987	enddo
2988	elseif(f2) then
2989	call make_node_layout(r) !!! bug (l) was wrong
2990	endif
2991	RETURN
2992	ELSE
2993	GOTO 1001
2994	ENDIF
2995	ENDIF
2996
2997
2998	! limit(1)=limit0(1)
2999	! limit(2)=limit0(2)
3000
3001	! CALL RING_L(R,doneit)
3002
3003	END SUBROUTINE THIN_LENS_resplit
3004
3005
3006
3007	SUBROUTINE RECUT_KIND7_one(C,lmax0,drift,ido,idc) ! A re-splitting routine
3008	IMPLICIT NONE
3009	TYPE(layout),POINTER :: L
3010	real(dp) lmax0
3011	TYPE(FIBRE),target :: C
3012	INTEGER I,f0,ido,idc
3013	logical(lp) drift,doit
3014
3015	if(associated(c%parent_layout)) then
3016	if(associated(c%parent_layout%parent_universe)) then
3017	l=>c%parent_layout%parent_universe%start
3018	do i=1,c%parent_layout%parent_universe%n
3019	call kill(l%t)
3020	l=>l%next
3021	enddo
3022	else
3023	call kill(c%parent_layout%t)
3024	endif
3025	endif
3026	if(drift.and.C%MAG%KIND==KIND1) then
3027	f0=nint(C%MAG%l/lmax0)
3028	if(f0==0) f0=1
3029	C%MAG%p%nst=f0
3030	C%MAGp%p%nst=C%MAG%p%nst
3031	call COPY(C%MAG,C%MAGP)
3032	endif
3033	if(mod(C%MAG%p%method,2)==0) then !!!
3034	doit=C%MAG%KIND==KIND7.or.(C%MAG%KIND==KIND2.and.C%MAG%p%method==2)
3035	if(associated(C%MAG%K16)) then
3036	doit=(C%MAG%K16%DRIFTKICK.and.C%MAG%p%method==2)
3037	endif
3038	if(associated(C%MAG%TP10)) then
3039	doit=(C%MAG%TP10%DRIFTKICK.and.C%MAG%p%method==2)
3040	endif
3041	IF(doit) then
3042	ido=ido+1
3043	! f0=nint(C%MAG%l/lmax0)
3044	f0=nint(C%MAG%l/lmax0/C%MAG%p%nst/2)
3045	if(C%MAG%p%method==6) f0=nint(C%MAG%l/lmax0/C%MAG%p%nst/4)
3046	if(f0==0) f0=1
3047	if(C%MAG%p%method==2) then
3048	C%MAG%p%nst=C%MAG%p%nstf02
3049	C%MAGp%p%nst=C%MAG%p%nst
3050	C%MAG%p%method=1
3051	C%MAGp%p%method=1
3052	elseif(C%MAG%p%method==4) then
3053	C%MAG%p%nst=C%MAG%p%nstf02
3054	C%MAGp%p%nst=C%MAG%p%nst
3055	C%MAG%p%method=3
3056	C%MAGp%p%method=3
3057	elseif(C%MAG%p%method==6) then
3058	C%MAG%p%nst=C%MAG%p%nstf04
3059	C%MAGp%p%nst=C%MAG%p%nst
3060	C%MAG%p%method=5
3061	C%MAGp%p%method=5
3062	endif
3063	call add(C%MAG,C%MAG%P%nmul,1,0.0_dp)
3064	call COPY(C%MAG,C%MAGP)
3065	if(C%MAG%KIND==KIND7) then
3066	C%MAG%t7%f=f0
3067	C%MAGp%t7%f=f0
3068	elseif(associated(C%MAG%K16)) then
3069	C%MAG%K16%f=f0
3070	C%MAGp%K16%f=f0
3071	elseif(associated(C%MAG%TP10)) then
3072	C%MAG%TP10%f=f0
3073	C%MAGp%TP10%f=f0
3074	else
3075	C%MAG%k2%f=f0
3076	C%MAGp%k2%f=f0
3077	endif
3078	ENDIF
3079	else !!!
3080	idc=idc+1
3081	f0=nint(C%MAG%l/lmax0/C%MAG%p%nst)
3082	if(f0>=1) then
3083	C%MAG%p%nst=C%MAG%p%nst*f0
3084	C%MAGp%p%nst=C%MAG%p%nst
3085	call add(C%MAG,C%MAG%P%nmul,1,0.0_dp)
3086	call COPY(C%MAG,C%MAGP)
3087	if(C%MAG%KIND==KIND7) then
3088	C%MAG%t7%f=f0*C%MAG%t7%f
3089	C%MAGp%t7%f=C%MAG%t7%f
3090	elseif(associated(C%MAG%K16)) then
3091	C%MAG%K16%f=f0*C%MAG%K16%f
3092	C%MAGp%K16%f=C%MAG%K16%f
3093	elseif(associated(C%MAG%TP10)) then
3094	C%MAG%TP10%f=f0*C%MAG%TP10%f
3095	C%MAGp%TP10%f=C%MAG%TP10%f
3096	else
3097	C%MAG%k2%f=f0*C%MAG%k2%f
3098	C%MAGp%k2%f=C%MAG%k2%f
3099	endif
3100	endif
3101	endif !!!
3102	END SUBROUTINE RECUT_KIND7_one
3103
3104
3105	SUBROUTINE RECUT_KIND7(R,lmax0,drift) ! A re-splitting routine
3106	IMPLICIT NONE
3107	TYPE(layout),target, intent(inout) :: R
3108	real(dp) lmax0
3109	TYPE(FIBRE),POINTER :: C
3110	INTEGER I,ido,idc
3111	logical(lp) drift
3112
3113	ido=0
3114	idc=0
3115
3116
3117	C=>R%START
3118	DO I=1,R%N
3119	call RECUT_KIND7_one(c,lmax0,drift,ido,idc)
3120	C=>C%NEXT
3121	ENDDO
3122	write(6,*) ido," elements changed to odd methods "
3123	write(6,*) idc," elements only "
3124	END SUBROUTINE RECUT_KIND7
3125
3126	SUBROUTINE ADD_SURVEY_INFO(R) ! A re-splitting routine
3127	IMPLICIT NONE
3128	TYPE(layout),target, intent(inout) :: R
3129	TYPE(FIBRE),POINTER :: C
3130	INTEGER I
3131	real(dp) b1
3132
3133
3134
3135	C=>R%START
3136	DO I=1,R%N
3137	if(C%mag%kind==kind3) then
3138	b1=C%mag%k3%thin_h_angle
3139	if(b1/=0.0_dp) then
3140	b1=b1/2.0_dp
3141	C%mag%k3%patch=my_true
3142	C%magp%k3%patch=my_true
3143	c%patch%patch=3
3144	c%patch%A_d=0.d0
3145	c%patch%B_d=0.d0
3146	c%patch%A_ANG=0.d0
3147	c%patch%B_ANG=0.d0
3148	c%patch%A_ANG(2)=b1
3149	c%patch%B_ANG(2)=b1
3150	endif
3151	endif
3152	C=>C%NEXT
3153	ENDDO
3154
3155	END SUBROUTINE ADD_SURVEY_INFO
3156
3157
3158
3159	SUBROUTINE check_bend(xl,ggi,rhoi,xbend1,gf,met) ! A re-splitting routine
3160	IMPLICIT NONE
3161	real(dp) xl,gg,ggi,rhoi,ar,ggb,co(7),xbend1,gf(7)
3162	integer i,met
3163
3164	gg=int(ggi)
3165	if(gg==0.d0) gg=1.d0
3166	co(3)=1.d0/6.0_dp
3167	co(5)= 0.2992989446749238e0_dp
3168	co(7)=0.2585213173527224e-1_dp
3169	ar=abs(rhoi)
3170	gf=0.0_dp
3171	! do i=3,7,2
3172	! ggb=((XL/gg) * ar)*(i)(XL/gg)*co(i) ! approximate residual orbit
3173	! if(xbend1<ggb) then
3174	! gf(i)=(ar*(i)co(i)/xbend1)*(one/(i+one))xl
3175	! write(6,*) i,gf(i)
3176	! pause
3177	! endif
3178	! enddo
3179
3180	do i=3,7,2
3181	ggb=((XL/gg) * ar)*(i-1)(XL/gg)co(i)twopi ! approximate residual orbit
3182	! if(xbend1<ggb) then
3183	gf(i)=(ar*(i-1)co(i)/xbend1/twopi)*(1.0_dp/(i+0.0_dp))xl
3184	! endif
3185	enddo
3186	gf(2)=gf(3)
3187	gf(4)=gf(5)*3
3188	gf(6)=gf(7)*7
3189
3190	met=2
3191
3192	if(gf(4)<gf(2)) met=4
3193	if(gf(6)<gf(4).and.gf(6)<gf(2)) met=6
3194	if(radiation_bend_split) met=2
3195	if(sixtrack_compatible) met=2
3196
3197	end SUBROUTINE check_bend
3198
3199	SUBROUTINE THIN_LENS_restart(R,fib,useknob) ! A re-splitting routine
3200	IMPLICIT NONE
3201	INTEGER NTE
3202	TYPE(layout),target, intent(inout) :: R
3203	INTEGER M1,M2,M3, MK1,MK2,MK3,nst_tot,ii !,limit0(2)
3204	type(fibre), OPTIONAL, target :: fib
3205	logical(lp), OPTIONAL :: useknob
3206	logical(lp) doit
3207	TYPE (fibre), POINTER :: C
3208	TYPE (layout), POINTER :: l
3209	! logical(lp) doneit
3210	nullify(C)
3211
3212	! CALL LINE_L(R,doneit)
3213
3214	if(associated(r%parent_universe)) then
3215	l=>r%parent_universe%start
3216	do ii=1,r%parent_universe%n
3217	call kill(l%t)
3218	l=>l%next
3219	enddo
3220	else
3221	call kill(r%t)
3222	endif
3223
3224
3225
3226
3227
3228	M1=0
3229	M2=0
3230	M3=0
3231	MK1=0
3232	MK2=0
3233	MK3=0
3234
3235
3236	r%NTHIN=0
3237
3238	nst_tot=0
3239	C=>R%START
3240	do ii=1,r%n ! WHILE(ASSOCIATED(C))
3241	doit=.true.
3242	if(present(useknob)) then
3243	if(useknob) then
3244	doit=c%magp%knob
3245	else
3246	doit=.not.c%magp%knob
3247	endif
3248	endif
3249	if(present(fib)) then
3250	doit=associated(c,fib)
3251	endif
3252	if(.not.doit) then
3253	NST_tot=NST_tot+C%MAG%P%nst
3254	C=>C%NEXT
3255	cycle
3256	endif
3257
3258
3259	doit=(C%MAG%KIND==kind1.and.C%MAG%KIND==kind2.or.C%MAG%KIND==kind5.or.C%MAG%KIND==kind4)
3260	doit=DOIT.OR.(C%MAG%KIND==kind6.or.C%MAG%KIND==kind7)
3261	DOIT=DOIT.OR.(C%MAG%KIND==kind10.or.C%MAG%KIND==kind16)
3262	DOIT=DOIT.OR.(C%MAG%KIND==kind17)
3263
3264
3265	if(doit) then
3266
3267	M1=M1+1
3268	NTE=1
3269	C%MAG%P%NST=NTE
3270	MK1=MK1+NTE
3271	C%MAG%P%METHOD=2
3272
3273	r%NTHIN=r%NTHIN+1 !C%MAG%NST
3274
3275	call add(C%MAG,C%MAG%P%nmul,1,0.0_dp)
3276	call COPY(C%MAG,C%MAGP)
3277	else
3278	if(C%MAG%KIND/=kindpa) then
3279	C%MAG%P%NST=1
3280	if(associated(C%MAG%bn))call add(C%MAG,C%MAG%P%nmul,1,0.0_dp)
3281	call COPY(C%MAG,C%MAGP)
3282	endif
3283	endif
3284
3285	NST_tot=NST_tot+C%MAG%P%nst
3286	C=>C%NEXT
3287	enddo
3288
3289
3290	write(6,*) "Present of cutable Elements ",r%NTHIN
3291	write(6,*) "METHOD 2 ",M1,MK1
3292	write(6,*) "METHOD 4 ",M2,MK2
3293	write(6,*) "METHOD 6 ",M3,MK3
3294	write(6,*) "number of Slices ", MK1+MK2+MK3
3295	write(6,*) "Total NST ", NST_tot
3296
3297
3298
3299	! CALL RING_L(R,doneit)
3300
3301	END SUBROUTINE THIN_LENS_restart
3302
3303
3304	SUBROUTINE print_bn_an(r,n,title,filename)
3305	implicit none
3306	type(layout),target,intent(inout) ::r
3307	character(*) filename
3308	type(fibre),pointer ::p
3309	integer n,i,mf,j,ntot
3310	character(*) title
3311
3312	ntot=0
3313	call kanalnummer(mf)
3314	open(unit=mf,file=filename)
3315	p=>r%start
3316	write(mf,'(a120)') title
3317	write(mf,*) n
3318	do i=1,r%n
3319
3320	if(associated(p%mag%an)) then
3321	ntot=ntot+1
3322	write(mf,*) min(n,p%mag%p%nmul),p%mag%name
3323	do j=1,min(n,p%mag%p%nmul)
3324	write(mf,*)j,p%mag%bn(j),p%mag%an(j)
3325	enddo
3326	endif
3327	p=>p%next
3328	enddo
3329
3330
3331	close(mf)
3332
3333	write(6,*) ntot," magnets settings saved to maximum order ",n
3334
3335	end SUBROUTINE print_bn_an
3336
3337	SUBROUTINE read_bn_an(r,filename)
3338	implicit none
3339	type(layout),target,intent(inout) ::r
3340	character(*) filename
3341	type(fibre),pointer ::p
3342	integer n,i,mf,j,nt,jt,ntot
3343	character(nlp) nom
3344	character*120 title
3345	real(dp), allocatable :: an(:),bn(:)
3346
3347	ntot=0
3348	call kanalnummer(mf)
3349	open(unit=mf,file=filename)
3350
3351	p=>r%start
3352	read(mf,'(a120)') title
3353	write(6,'(a120)') title
3354	read(mf,*) n
3355	allocate(an(n),bn(n))
3356	an=0.0_dp;bn=0.0_dp;
3357
3358	do i=1,r%n
3359
3360	if(associated(p%mag%an)) then
3361	read(mf,*) nt,nom
3362	call context(nom)
3363
3364	do j=1,nt
3365	read(mf,*)jt,bn(j),an(j)
3366	enddo
3367
3368	ntot=ntot+1
3369	do j=nt,1,-1
3370	call ADD(p,j,0,bn(j))
3371	call ADD(p,-j,0,an(j))
3372	enddo
3373	endif ! associated
3374	p=>p%next
3375	enddo
3376
3377	write(6,*) ntot," magnets settings read"
3378
3379	close(mf)
3380	deallocate(an,bn)
3381
3382	end SUBROUTINE read_bn_an
3383
3384	! THIN LENS EXAMPLE
3385
3386	SUBROUTINE assign_one_aperture(L,pos,kindaper,R,X,Y,dx,dy)
3387	IMPLICIT NONE
3388	TYPE(LAYOUT),TARGET :: L
3389	integer pos,kindaper
3390	REAL(DP) R(:),X,Y,dx,dy
3391	type(fibre), pointer :: P
3392
3393	call move_to(L,p,pos)
3394
3395	if(.NOT.ASSOCIATED(P%MAG%p%aperture)) THEN
3396	call alloc(P%MAG%p%aperture)
3397	call alloc(P%MAGP%p%aperture)
3398	ENDIF
3399	if(kindaper/=0) then
3400	P%MAG%p%aperture%kind = kindaper
3401	P%MAGP%p%aperture%kind = kindaper
3402	P%MAG%p%aperture%r = R
3403	P%MAG%p%aperture%x = X
3404	P%MAG%p%aperture%y = y
3405	P%MAG%p%aperture%dx = dX
3406	P%MAG%p%aperture%dy = dy
3407	P%MAGP%p%aperture%r = R
3408	P%MAGP%p%aperture%x = X
3409	P%MAGP%p%aperture%y = y
3410	P%MAGP%p%aperture%dx = dX
3411	P%MAGP%p%aperture%dy = dy
3412	endif
3413
3414	end SUBROUTINE assign_one_aperture
3415
3416	SUBROUTINE TURN_OFF_ONE_aperture(R,pos)
3417	IMPLICIT NONE
3418	TYPE(LAYOUT),TARGET :: R
3419	integer pos
3420	type(fibre), pointer :: P
3421
3422	call move_to(r,p,pos)
3423
3424	if(ASSOCIATED(P%MAG%p%aperture)) THEN
3425	P%MAG%p%aperture%kind = -P%MAG%p%aperture%kind
3426	P%MAGP%p%aperture%kind = P%MAG%p%aperture%kind
3427	ENDIF
3428
3429	end SUBROUTINE TURN_OFF_ONE_aperture
3430
3431	SUBROUTINE REVERSE_BEAM_LINE(R, changeanbn )
3432	IMPLICIT NONE
3433	TYPE(LAYOUT),TARGET :: R
3434	integer J,I
3435	type(fibre), pointer :: P
3436	logical(lp), optional:: changeanbn
3437	logical(lp) changeanbn0
3438
3439	changeanbn0=my_true
3440	if(present(changeanbn)) changeanbn0=changeanbn
3441
3442	p=>r%start
3443	do i=1,r%n
3444	p%dir=-1
3445	if(changeanbn0) then
3446	if(associated(p%mag%an)) then
3447	do j=1,p%mag%p%nmul
3448	p%mag%bn(j)=-p%magp%bn(j)
3449	p%mag%an(j)=-p%magp%an(j)
3450	p%magp%bn(j)=-p%magp%bn(j)
3451	p%magp%an(j)=-p%magp%an(j)
3452	enddo
3453	if(abs(abs(p%mag%bn(1))-abs(p%mag%p%b0)).gt.1e-11_dp.or. &
3454	abs(p%mag%p%b0).lt.1e-11_dp) then
3455	p%mag%bn(1)=-p%magp%bn(1)
3456	p%mag%an(1)=-p%magp%an(1)
3457	p%magp%bn(1)=-p%magp%bn(1)
3458	p%magp%an(1)=-p%magp%an(1)
3459	endif
3460	if(p%mag%p%nmul>0) call add(p,1,1,0.0_dp)
3461	endif
3462	if(associated(p%mag%volt)) p%mag%volt=-p%mag%volt
3463	if(associated(p%magp%volt)) p%magp%volt=-p%magp%volt
3464	endif
3465	P=>P%next
3466	ENDDO
3467	end SUBROUTINE REVERSE_BEAM_LINE
3468
3469	SUBROUTINE PUTFRINGE(R, changeanbn )
3470	IMPLICIT NONE
3471	TYPE(LAYOUT),TARGET :: R
3472	integer I
3473	type(fibre), pointer :: P
3474	logical(lp) changeanbn
3475
3476	p=>r%start
3477	do i=1,r%n
3478	p%mag%p%PERMFRINGE =changeanbn
3479	p%magP%p%PERMFRINGE=changeanbn
3480	P=>P%next
3481	ENDDO
3482	end SUBROUTINE PUTFRINGE
3483
3484	SUBROUTINE PUTbend_fringe(R, changeanbn )
3485	IMPLICIT NONE
3486	TYPE(LAYOUT),TARGET :: R
3487	integer I
3488	type(fibre), pointer :: P
3489	logical(lp) changeanbn
3490
3491	p=>r%start
3492	do i=1,r%n
3493	if(p%mag%p%b0/=0.0_dp) then
3494	p%mag%p%bend_fringe =changeanbn
3495	p%magp%p%bend_fringe =changeanbn
3496	write(6,*) P%mag%name, " changed to ",changeanbn
3497	endif
3498	P=>P%next
3499	ENDDO
3500	end SUBROUTINE PUTbend_fringe
3501
3502
3503	SUBROUTINE MESS_UP_ALIGNMENT(R,SIG,cut)
3504	use gauss_dis
3505	IMPLICIT NONE
3506	TYPE(LAYOUT),TARGET :: R
3507	integer J,I
3508	type(fibre), pointer :: P
3509	REAL(DP) SIG(:),X,MIS(6),cut
3510
3511
3512	p=>r%start
3513	do i=1,r%n
3514
3515	IF(P%MAG%KIND/=KIND0.AND.P%MAG%KIND/=KIND1) THEN
3516	DO J=1,6
3517	call GRNF(X,cut)
3518	MIS(J)=X*SIG(J)
3519	ENDDO
3520	call MISALIGN_FIBRE(p,mis)
3521	ENDIF
3522	P=>P%NEXT
3523	ENDDO
3524	end SUBROUTINE MESS_UP_ALIGNMENT
3525
3526
3527	! CALL MESS_UP_ALIGNMENT_name(my_ring,name,i1,i2,SIG,cut)
3528
3529	subroutine MESS_UP_ALIGNMENT_name(R,nom,iseed,full,sig,cut,pr)
3530	use gauss_dis
3531	IMPLICIT NONE
3532	TYPE(layout),target, intent(inout):: R
3533	integer iseed,j,ic,i
3534	character(nlp) nom
3535	type(fibre), pointer :: p
3536	logical(lp) f1,full,pr
3537	real(dp) cut,sig(6),mis(6),x,taxi(6)
3538
3539	if(iseed/=0) call gaussian_seed(iseed)
3540
3541	taxi=0.d0
3542	call context(nom)
3543
3544	ic=0
3545
3546	p=>r%start
3547	do i=1,r%n
3548	f1=.false.
3549
3550	if(full) then
3551	f1=(p%mag%name ==nom )
3552	else
3553	f1=( index(p%mag%name,nom(1:len_trim(nom))) >0)
3554	endif
3555
3556
3557	if(f1) then
3558	ic=ic+1
3559	DO J=1,6
3560	call GRNF(X,cut)
3561	MIS(J)=X*SIG(J)
3562	taxi(j)=taxi(j)+abs(MIS(J))
3563	ENDDO
3564	call MISALIGN_FIBRE(p,mis)
3565	endif
3566
3567	if(f1) then
3568	ic=ic+1
3569	if(pr) write(6,*) p%mag%name
3570	endif
3571
3572	p=>P%next
3573	enddo
3574
3575	if(pr) then
3576	write(6,*) ic," Magnets modified "
3577
3578
3579
3580	write(6,*) ic," Magnets misaligned "
3581	taxi=taxi/ic
3582
3583	write(6,'(a21,3(1x,E15.8))') " <\|displacements\|> = ",taxi(1:3)
3584	write(6,'(a21,3(1x,E15.8))') " <\|rotations\|> = ",taxi(4:6)
3585	endif
3586
3587	end subroutine MESS_UP_ALIGNMENT_name
3588
3589	subroutine Sigma_of_alignment(r,sig,ave)
3590	IMPLICIT NONE
3591	TYPE(layout),target, intent(inout):: R
3592	integer i,j,is(6)
3593	type(fibre), pointer :: p
3594	real(dp) sig(6),ave(6)
3595	character*23 lab(6)
3596	lab(1)=" Dx Average and Sigma "
3597	lab(2)=" Dy Average and Sigma "
3598	lab(3)=" Dz Average and Sigma "
3599	lab(4)=" Dax Average and Sigma "
3600	lab(5)=" Day Average and Sigma "
3601	lab(6)=" Daz Average and Sigma "
3602
3603	AVE=0.0_dp
3604	SIG=0.0_dp
3605	p=>r%start
3606	is=0
3607	do i=1,r%n
3608	do j=1,3
3609	if(p%chart%D_IN(j)/=0.0_dp) then
3610	is(j)=is(j)+1
3611	ave(j)=p%chart%D_IN(j)+ave(j)
3612	sig(j)=p%chart%D_IN(j)**2+sig(j)
3613	endif
3614	enddo
3615	do j=4,6
3616	if(p%chart%ANG_IN(j)/=0.0_dp) then
3617	is(j)=is(j)+1
3618	ave(j)=p%chart%ANG_IN(j)+ave(j)
3619	sig(j)=p%chart%ANG_IN(j)**2+sig(j)
3620	endif
3621	enddo
3622	p=>p%next
3623	enddo
3624
3625	do i=1,6
3626	if(is(i)/=0) then
3627	ave(i)=ave(i)/is(i)
3628	sig(i)=sig(i)/is(i)
3629	sig(i)=sqrt(sig(i)-ave(i)**2)
3630	write(6,*) is(i), " Magnets misaligned "
3631	write(6,"(1x,a23,2(1x,E15.8))") lab(i),ave(i),sig(i)
3632	endif
3633	enddo
3634
3635	end subroutine Sigma_of_alignment
3636
3637
3638
3639
3640	SUBROUTINE dyn_aper(L,x_in,n_in,ang_in,ang_out,del_in,dlam,pos,nturn,ite,state,mf)
3641	IMPLICIT NONE
3642	type(layout),target, intent(inout) :: L
3643	real(dp) x(6)
3644	REAL(DP) x_in,del_in,closed(6),r(6),rt(6)
3645	REAL(DP) lamT,lams,lamu,dlam,DLAMT,DX,ang,ang_in,ang_out
3646	integer pos,nturn,i,st,ite,ic,mf,n_in,j_in
3647	TYPE(INTERNAL_STATE) STATE
3648	!
3649	! TYPE(REAL_8) Y(6)
3650	! TYPE(DAMAP) ID
3651	! TYPE(NORMALFORM) NORM
3652
3653	closed=0.0_dp
3654	! STATE=STATE+NOCAVITY0
3655	if(state%nocavity) closed(5)=del_in
3656
3657	CALL FIND_ORBIT(L,CLOSED,pos,STATE,1e-5_dp)
3658	write(6,*) "closed orbit "
3659	write(6,*) CLOSED
3660	write(mf,201) closed
3661	ang= (ang_out-ang_in)/n_in
3662	lamt=1.0_dp
3663	do j_in=0,n_in
3664
3665	x=0.0_dp
3666	x(1)=x_incos(j_inang+ang_in)
3667	x(3)=x_insin(j_inang+ang_in)
3668	x(5)=del_in
3669
3670
3671	dx=0.3_dp
3672
3673	r=0.0_dp;rt=0.0_dp;
3674	lams=0.0_dp
3675	lamu=0.0_dp
3676
3677	DLAMT=DX
3678
3679	! lamt=ONE
3680	ic=0
3681	do while(DLAMT>dlam.and.ic<ite)
3682
3683	ic=ic+1
3684	R=0.0_dp;
3685	r(1:4)=lamt*x(1:4)
3686	if(state%nocavity) then
3687	rt=r+closed
3688	else
3689	rt=r+closed
3690	rt(5)=rt(5)+x(5)
3691	endif
3692
3693
3694	do i=1,nturn
3695	st=track_flag(L,rt,pos,state)
3696	if(st/=0) exit
3697	enddo
3698
3699	if(st/=0) then
3700	lamu=lamt
3701	lamt=(lams+lamt)/2.0_dp
3702	else
3703	lams=lamt
3704	IF(LAMU<DX) THEN
3705	lamt=DX+lamt
3706	ELSE
3707	lamt=(lamu+lamt)/2.0_dp
3708	ENDIF
3709	endif
3710	DLAMT=sqrt(x(1)2+x(3)2)*ABS(LAMU-LAMS)
3711	enddo
3712	write(6,) ic,(j_inang+ang_in)/twopi,lamSx(1),lamSx(3)
3713
3714	write(mf,202) lamSx(1),lamSx(3),lamSx(1)+closed(1),lamSx(3)+closed(3),DLAMT
3715	lamt=lamt*0.8_dp
3716	enddo
3717	201 FORMAT(6(1X,D18.11))
3718	202 FORMAT(5(1X,D18.11))
3719
3720	end SUBROUTINE dyn_aper
3721
3722
3723
3724	SUBROUTINE dyn_aperalex(L,x_in,del_in,dx,dlam,pos,nturn,ite,state,mf,targ_tune,fixp)
3725	IMPLICIT NONE
3726	type(layout),target, intent(inout) :: L
3727	real(dp) x(6)
3728	REAL(DP) x_in,del_in,closed(6),r(6),rt(6),targ_tune(2)
3729	REAL(DP) lamT,lams,lamu,dlam,DLAMT,DX,ang,ang_in,ang_out
3730	integer pos,nturn,i,st,ite,ic,mf,n_in,j_in
3731	TYPE(INTERNAL_STATE) STATE
3732	logical(lp) fixp
3733	!
3734	! TYPE(REAL_8) Y(6)
3735	! TYPE(DAMAP) ID
3736	! TYPE(NORMALFORM) NORM
3737
3738	closed=0.0_dp
3739	! STATE=STATE+NOCAVITY0
3740
3741	! if(state%nocavity)
3742	closed(5)=del_in
3743	if(fixp) then
3744	CALL FIND_ORBIT(L,CLOSED,pos,STATE,1e-5_dp)
3745	write(6,*) "closed orbit "
3746	write(6,*) CLOSED
3747	else
3748	closed(1:4)=0.0_dp
3749	closed(6)=0.0_dp
3750	endif
3751	! write(mf,201) closed
3752	n_in=1
3753	ang_in=pi/4.0_dp
3754	ang_out=pi/4.0_dp
3755	ang= (ang_out-ang_in)/n_in
3756	lamt=1.0_dp
3757	j_in=0
3758	! do j_in=0,n_in
3759
3760	x=0.0_dp
3761	x(1)=x_incos(j_inang+ang_in)
3762	x(3)=x_insin(j_inang+ang_in)
3763	! x(5)=del_in
3764
3765
3766
3767	r=0.0_dp;rt=0.0_dp;
3768	lams=0.0_dp
3769	lamu=0.0_dp
3770
3771	DLAMT=DX
3772
3773	! lamt=ONE
3774	ic=0
3775	do while(DLAMT>dlam.and.ic<ite)
3776
3777	ic=ic+1
3778	R=0.0_dp;
3779	r(1:4)=lamt*x(1:4)
3780	if(state%nocavity) then
3781	if(fixp) then
3782	rt=r+closed
3783	else
3784	rt(1:4)=r(1:4)
3785	rt(6)=0.0_dp
3786	rt(5)=del_in
3787	endif
3788	else
3789	if(fixp) then
3790	rt=r+closed
3791	rt(5)=rt(5)+del_in
3792	else
3793	rt=r
3794	rt(5)=rt(5)+del_in
3795	endif
3796	endif
3797
3798	! write(6,*) rt
3799	do i=1,nturn
3800	st=track_flag(L,rt,pos,state)
3801	if(st/=0) exit
3802	enddo
3803	! write(6,*) i,check_stable
3804	! write(6,*) rt
3805	! pause
3806	if(st/=0) then
3807	lamu=lamt
3808	lamt=(lams+lamt)/2.0_dp
3809	! write(mf,*) "unstable ",lamt
3810
3811	else
3812	lams=lamt
3813	IF(LAMU<DX) THEN
3814	lamt=DX+lamt
3815	ELSE
3816	lamt=(lamu+lamt)/2.0_dp
3817	ENDIF
3818	! write(mf,*) "stable ",lamt
3819
3820	endif
3821	DLAMT=sqrt(x(1)2+x(3)2)*ABS(LAMU-LAMS)
3822	! write(mf,) "dlamt ",dlamt,sqrt(x(1)2+x(3)*2)
3823
3824	enddo
3825	write(6,) ic,lamSx(1) !,lamSx(3),(j_inang+ang_in)/twopi
3826
3827	write(mf,202) targ_tune,lamSx(1) !,lamSx(3),lamSx(1)+closed(1),lamSx(3)+closed(3),DLAMT
3828	lamt=lamt*0.8_dp
3829	! enddo
3830	201 FORMAT(6(1X,D18.11))
3831	202 FORMAT(3(1X,D18.11))
3832
3833	end SUBROUTINE dyn_aperalex
3834
3835
3836
3837
3838	end module S_fitting

Note: See TracBrowser for help on using the repository browser.

Context Navigation

source: PSPA/madxPSPA/libs/ptc/src/So_fitting.f90 @ 430

Download in other formats: