2015-03-04 17:07:15 +00:00
										 
									 
								 
							 | 
							
								
									
										
									
								
							 | 
							
								
							 | 
							
							
								subroutine afc9(c3a,npts,fsample,a,syncpk)
							 | 
						
					
						
							
								
									
										
										
										
											2013-08-10 15:29:55 +00:00
										 
									 
								 
							 | 
							
								
									
										
									
								
							 | 
							
								
							 | 
							
							
								
							 | 
						
					
						
							
								
									
										
										
										
											2015-04-22 17:48:03 +00:00
										 
									 
								 
							 | 
							
								
									
										
									
								
							 | 
							
								
							 | 
							
							
								  parameter (NZ2=1512)
							 | 
						
					
						
							| 
								
							 | 
							
								
							 | 
							
								
							 | 
							
							
								  complex c3a(0:NZ2-1)
							 | 
						
					
						
							| 
								
							 | 
							
								
							 | 
							
								
							 | 
							
							
								  complex c3(0:NZ2-1)
							 | 
						
					
						
							
								
									
										
										
										
											2013-08-10 15:29:55 +00:00
										 
									 
								 
							 | 
							
								
									
										
									
								
							 | 
							
								
							 | 
							
							
								  real a(3),deltaa(3)
							 | 
						
					
						
							| 
								
							 | 
							
								
							 | 
							
								
							 | 
							
							
								
							 | 
						
					
						
							| 
								
							 | 
							
								
							 | 
							
								
							 | 
							
							
								  a(1)=0.                                   !f0
							 | 
						
					
						
							| 
								
							 | 
							
								
							 | 
							
								
							 | 
							
							
								  a(2)=0.                                   !f1
							 | 
						
					
						
							| 
								
							 | 
							
								
							 | 
							
								
							 | 
							
							
								  a(3)=0.                                   !f2
							 | 
						
					
						
							
								
									
										
										
										
											2015-03-04 17:07:15 +00:00
										 
									 
								 
							 | 
							
								
									
										
									
								
							 | 
							
								
							 | 
							
							
								  deltaa(1)=1.736
							 | 
						
					
						
							| 
								
							 | 
							
								
							 | 
							
								
							 | 
							
							
								  deltaa(2)=1.736
							 | 
						
					
						
							| 
								
							 | 
							
								
							 | 
							
								
							 | 
							
							
								  deltaa(3)=1.0
							 | 
						
					
						
							
								
									
										
										
										
											2013-08-10 15:29:55 +00:00
										 
									 
								 
							 | 
							
								
									
										
									
								
							 | 
							
								
							 | 
							
							
								  nterms=3
							 | 
						
					
						
							| 
								
							 | 
							
								
							 | 
							
								
							 | 
							
							
								
							 | 
						
					
						
							| 
								
							 | 
							
								
							 | 
							
								
							 | 
							
							
								! Start the iteration
							 | 
						
					
						
							| 
								
							 | 
							
								
							 | 
							
								
							 | 
							
							
								  chisqr=0.
							 | 
						
					
						
							| 
								
							 | 
							
								
							 | 
							
								
							 | 
							
							
								  chisqr0=1.e6
							 | 
						
					
						
							
								
									
										
										
										
											2015-03-04 17:07:15 +00:00
										 
									 
								 
							 | 
							
								
									
										
									
								
							 | 
							
								
							 | 
							
							
								  c3=c3a
							 | 
						
					
						
							| 
								
							 | 
							
								
							 | 
							
								
							 | 
							
							
								  a3=a(3)
							 | 
						
					
						
							| 
								
							 | 
							
								
							 | 
							
								
							 | 
							
							
								  do iter=1,4
							 | 
						
					
						
							
								
									
										
										
										
											2013-08-10 15:29:55 +00:00
										 
									 
								 
							 | 
							
								
									
										
									
								
							 | 
							
								
							 | 
							
							
								     do j=1,nterms
							 | 
						
					
						
							
								
									
										
										
										
											2015-03-05 17:13:00 +00:00
										 
									 
								 
							 | 
							
								
									
										
									
								
							 | 
							
								
							 | 
							
							
								        if(a(3).ne.a3) call shft(c3a,a(3),a3,c3)
							 | 
						
					
						
							
								
									
										
										
										
											2013-08-10 15:29:55 +00:00
										 
									 
								 
							 | 
							
								
									
										
									
								
							 | 
							
								
							 | 
							
							
								        chisq1=fchisq(c3,npts,fsample,a)
							 | 
						
					
						
							| 
								
							 | 
							
								
							 | 
							
								
							 | 
							
							
								        fn=0.
							 | 
						
					
						
							| 
								
							 | 
							
								
							 | 
							
								
							 | 
							
							
								        delta=deltaa(j)
							 | 
						
					
						
							| 
								
							 | 
							
								
							 | 
							
								
							 | 
							
							
								10      a(j)=a(j)+delta
							 | 
						
					
						
							
								
									
										
										
										
											2015-03-05 17:13:00 +00:00
										 
									 
								 
							 | 
							
								
									
										
									
								
							 | 
							
								
							 | 
							
							
								        if(a(3).ne.a3) call shft(c3a,a(3),a3,c3)
							 | 
						
					
						
							
								
									
										
										
										
											2013-08-10 15:29:55 +00:00
										 
									 
								 
							 | 
							
								
									
										
									
								
							 | 
							
								
							 | 
							
							
								        chisq2=fchisq(c3,npts,fsample,a)
							 | 
						
					
						
							| 
								
							 | 
							
								
							 | 
							
								
							 | 
							
							
								        if(chisq2.eq.chisq1) go to 10
							 | 
						
					
						
							| 
								
							 | 
							
								
							 | 
							
								
							 | 
							
							
								        if(chisq2.gt.chisq1) then
							 | 
						
					
						
							| 
								
							 | 
							
								
							 | 
							
								
							 | 
							
							
								           delta=-delta                      !Reverse direction
							 | 
						
					
						
							| 
								
							 | 
							
								
							 | 
							
								
							 | 
							
							
								           a(j)=a(j)+delta
							 | 
						
					
						
							| 
								
							 | 
							
								
							 | 
							
								
							 | 
							
							
								           tmp=chisq1
							 | 
						
					
						
							| 
								
							 | 
							
								
							 | 
							
								
							 | 
							
							
								           chisq1=chisq2
							 | 
						
					
						
							| 
								
							 | 
							
								
							 | 
							
								
							 | 
							
							
								           chisq2=tmp
							 | 
						
					
						
							| 
								
							 | 
							
								
							 | 
							
								
							 | 
							
							
								        endif
							 | 
						
					
						
							| 
								
							 | 
							
								
							 | 
							
								
							 | 
							
							
								20      fn=fn+1.0
							 | 
						
					
						
							| 
								
							 | 
							
								
							 | 
							
								
							 | 
							
							
								        a(j)=a(j)+delta
							 | 
						
					
						
							
								
									
										
										
										
											2015-03-05 17:13:00 +00:00
										 
									 
								 
							 | 
							
								
									
										
									
								
							 | 
							
								
							 | 
							
							
								        if(a(3).ne.a3) call shft(c3a,a(3),a3,c3)
							 | 
						
					
						
							
								
									
										
										
										
											2013-08-10 15:29:55 +00:00
										 
									 
								 
							 | 
							
								
									
										
									
								
							 | 
							
								
							 | 
							
							
								        chisq3=fchisq(c3,npts,fsample,a)
							 | 
						
					
						
							| 
								
							 | 
							
								
							 | 
							
								
							 | 
							
							
								        if(chisq3.lt.chisq2) then
							 | 
						
					
						
							| 
								
							 | 
							
								
							 | 
							
								
							 | 
							
							
								           chisq1=chisq2
							 | 
						
					
						
							| 
								
							 | 
							
								
							 | 
							
								
							 | 
							
							
								           chisq2=chisq3
							 | 
						
					
						
							| 
								
							 | 
							
								
							 | 
							
								
							 | 
							
							
								           go to 20
							 | 
						
					
						
							| 
								
							 | 
							
								
							 | 
							
								
							 | 
							
							
								        endif
							 | 
						
					
						
							| 
								
							 | 
							
								
							 | 
							
								
							 | 
							
							
								
							 | 
						
					
						
							| 
								
							 | 
							
								
							 | 
							
								
							 | 
							
							
								! Find minimum of parabola defined by last three points
							 | 
						
					
						
							| 
								
							 | 
							
								
							 | 
							
								
							 | 
							
							
								        delta=delta*(1./(1.+(chisq1-chisq2)/(chisq3-chisq2))+0.5)
							 | 
						
					
						
							| 
								
							 | 
							
								
							 | 
							
								
							 | 
							
							
								        a(j)=a(j)-delta
							 | 
						
					
						
							
								
									
										
										
										
											2015-03-04 17:07:15 +00:00
										 
									 
								 
							 | 
							
								
									
										
									
								
							 | 
							
								
							 | 
							
							
								        if(j.lt.3) deltaa(j)=deltaa(j)*fn/3.
							 | 
						
					
						
							
								
									
										
										
										
											2015-02-27 15:38:24 +00:00
										 
									 
								 
							 | 
							
								
									
										
									
								
							 | 
							
								
							 | 
							
							
								!        write(*,4000) iter,j,a,-chisq2
							 | 
						
					
						
							| 
								
							 | 
							
								
							 | 
							
								
							 | 
							
							
								!4000    format(i1,i2,3f10.4,f11.3)
							 | 
						
					
						
							
								
									
										
										
										
											2013-08-10 15:29:55 +00:00
										 
									 
								 
							 | 
							
								
									
										
									
								
							 | 
							
								
							 | 
							
							
								     enddo
							 | 
						
					
						
							
								
									
										
										
										
											2015-03-05 17:13:00 +00:00
										 
									 
								 
							 | 
							
								
									
										
									
								
							 | 
							
								
							 | 
							
							
								     if(a(3).ne.a3) call shft(c3a,a(3),a3,c3)
							 | 
						
					
						
							
								
									
										
										
										
											2013-08-10 15:29:55 +00:00
										 
									 
								 
							 | 
							
								
									
										
									
								
							 | 
							
								
							 | 
							
							
								     chisqr=fchisq(c3,npts,fsample,a)
							 | 
						
					
						
							
								
									
										
										
										
											2016-01-08 20:55:17 +00:00
										 
									 
								 
							 | 
							
								
									
										
									
								
							 | 
							
								
							 | 
							
							
								!     write(*,4000) 0,0,a,-chisqr
							 | 
						
					
						
							
								
									
										
										
										
											2015-02-27 15:38:24 +00:00
										 
									 
								 
							 | 
							
								
									
										
									
								
							 | 
							
								
							 | 
							
							
								     if(chisqr/chisqr0.gt.0.99) exit
							 | 
						
					
						
							
								
									
										
										
										
											2013-08-10 15:29:55 +00:00
										 
									 
								 
							 | 
							
								
									
										
									
								
							 | 
							
								
							 | 
							
							
								     chisqr0=chisqr
							 | 
						
					
						
							| 
								
							 | 
							
								
							 | 
							
								
							 | 
							
							
								  enddo
							 | 
						
					
						
							| 
								
							 | 
							
								
							 | 
							
								
							 | 
							
							
								
							 | 
						
					
						
							| 
								
							 | 
							
								
							 | 
							
								
							 | 
							
							
								  syncpk=-chisqr
							 | 
						
					
						
							
								
									
										
										
										
											2015-03-04 17:07:15 +00:00
										 
									 
								 
							 | 
							
								
									
										
									
								
							 | 
							
								
							 | 
							
							
								  c3a=c3
							 | 
						
					
						
							
								
									
										
										
										
											2016-01-08 20:55:17 +00:00
										 
									 
								 
							 | 
							
								
									
										
									
								
							 | 
							
								
							 | 
							
							
								!  write(*,4001) a,syncpk
							 | 
						
					
						
							
								
									
										
										
										
											2015-02-27 15:38:24 +00:00
										 
									 
								 
							 | 
							
								
									
										
									
								
							 | 
							
								
							 | 
							
							
								!4001 format(3x,3f10.4,f11.3)
							 | 
						
					
						
							
								
									
										
										
										
											2013-08-10 15:29:55 +00:00
										 
									 
								 
							 | 
							
								
									
										
									
								
							 | 
							
								
							 | 
							
							
								
							 | 
						
					
						
							| 
								
							 | 
							
								
							 | 
							
								
							 | 
							
							
								  return
							 | 
						
					
						
							| 
								
							 | 
							
								
							 | 
							
								
							 | 
							
							
								end subroutine afc9
							 | 
						
					
						
							
								
									
										
										
										
											2015-03-05 17:13:00 +00:00
										 
									 
								 
							 | 
							
								
									
										
									
								
							 | 
							
								
							 | 
							
							
								
							 | 
						
					
						
							| 
								
							 | 
							
								
							 | 
							
								
							 | 
							
							
								subroutine shft(c3a,a3a,a3,c3)
							 | 
						
					
						
							| 
								
							 | 
							
								
							 | 
							
								
							 | 
							
							
								  complex c3a(0:1359)
							 | 
						
					
						
							| 
								
							 | 
							
								
							 | 
							
								
							 | 
							
							
								  complex  c3(0:1359)
							 | 
						
					
						
							| 
								
							 | 
							
								
							 | 
							
								
							 | 
							
							
								
							 | 
						
					
						
							| 
								
							 | 
							
								
							 | 
							
								
							 | 
							
							
								  a3=a3a
							 | 
						
					
						
							| 
								
							 | 
							
								
							 | 
							
								
							 | 
							
							
								  n=nint(a3)
							 | 
						
					
						
							| 
								
							 | 
							
								
							 | 
							
								
							 | 
							
							
								  c3=cshift(c3a,n)
							 | 
						
					
						
							| 
								
							 | 
							
								
							 | 
							
								
							 | 
							
							
								  if(n.gt.0) c3(1360-n:)=0.0
							 | 
						
					
						
							| 
								
							 | 
							
								
							 | 
							
								
							 | 
							
							
								  if(n.lt.0) c3(:n-1)=0.0
							 | 
						
					
						
							| 
								
							 | 
							
								
							 | 
							
								
							 | 
							
							
								
							 | 
						
					
						
							| 
								
							 | 
							
								
							 | 
							
								
							 | 
							
							
								  return
							 | 
						
					
						
							| 
								
							 | 
							
								
							 | 
							
								
							 | 
							
							
								end subroutine shft
							 |