c #################################################################        
c program mca.f  
c Ajustement par moindres carres   
c    n   : number of involved observations        
c    l   : number of parameter to be estimated
c    k1max: total number of waves    
c ##################################################################         
       
	implicit none        
	integer n1,k1,k2
	parameter(n1=50000,k1=20,k2=2*k1+4) 	       
	double precision za(n1,k2),zp(n1),pp(n1),xr(n1), xs(k2),xdx(k2),cor(k2,k2)  
                  	
	double precision TAB(20)
	
        double precision DATG(n1),DAT1(n1),X(n1),Y(n1),EX(n1),EY(n1)
	
	double precision PERIOD(k1),A(K1),AP(K1),DAP(K1), PH(K1),PH1,DPH(K1),ARGLS(5)  
      
      	double precision error,pi,dpi,date_deb,date_fin,an,date1,date2,
     &         span,sspan,datem,year,
     &         amp,da,phase,dphase,dt,dx,dy,per,date,xx,yy,chi2,X_mod,Y_mod,sx,sy
	      
        double precision rms,rms_r  
	
	integer arg(k1,6),cold,colv1,colv2,cole1,cole2,poly,npol,fac
          
        integer i,j,k,l,m,n,isel,ipris,m1,m2,m3,m4,m5,m6,k1max,nes,nmax,iw,maxc
	double precision t0
	character*80 kom1,filearg, filein,c_TT(20) 
	character*80 c_filein,c_cold,c_colv1,c_colv2,c_cole1,c_cole2,c_date_deb,c_date_fin,c_t0,c_span,c_sspan,c_T,c_poly,c_npol
	                                               
       
	pi=4.d0*datan(1.d0)        
	dpi=2.d0*pi              	

        call getarg(1,c_filein) 
	
	if (Len_trim(c_filein).eq.0) then 
	write(6,*) '            <<< mca  >>>          '
	write(6,*) ' Ajustement de termes periodiques ou polynomiaux dans une serie temporelle 1D ou 2D'
        write(6,*) ' mca [filein][cold][colv1][colv2][cole1][cole2][date_deb][date_fin][t0][span][sspan][T][poly][npol]'       
	write(6,*) ' [filein] .... fichier entree '       
	write(6,*) ' [cold] ...... colonne date   '       
	write(6,*) ' [colv1] ..... colonne des valeurs x ' 
	write(6,*) ' [colv2] ..... colonne des valeurs y (0 pour serie 1D) -colv2 pour prendre -y' 
	write(6,*) ' [cole1] ..... colonne des incertitudes sur x (0 pour ponderation uniforme)' 
	write(6,*) ' [cole2] ..... colonne des incertitudes sur y (0 pour ponderation uniforme)' 
	write(6,*) ' [date_deb] .. date de debut'
	write(6,*) ' [date_fin] .. date finale'
	write(6,*) ' [t0] ........ Initial epoch'
	write(6,*) ' [span] ...... periode estimation : 0 pour [date_deb, date_fin]'
	write(6,*) ' [sspan] ..... periode de glissement : 0 par defaut'
        write(6,*) ' [T] ......... Chaine T1_T2_T3_.... des periodes T1,T2,T3 .... des harmoniques a ajuster'
	write(6,*) "               T = 0 : n'ajuster aucune periode" 	 
	write(6,*) ' [poly] ...... Ajustement polynomial : 0 (non) / 1 (oui)'
	write(6,*) ' [npol] ...... degre du polynonme a ajuster : 0 par defaut'
	write(6,*) ' Resultats dans mca.res & residus dans mca.resi'
	write(6,*) '  '
	write(6,*) ' exemple : mca extract.res 1 2 3 0 0 50000 55000 51544 0 0 365.25_433 1 1' 
        stop 
	endif
        
	filein=c_filein      	
	call getarg(2,c_cold)
	call getarg(3,c_colv1)
	call getarg(4,c_colv2)
	call getarg(5,c_cole1)
	call getarg(6,c_cole2)
	call getarg(7,c_date_deb)
	call getarg(8,c_date_fin)
	call getarg(9,c_t0)
	call getarg(10,c_span)
	call getarg(11,c_sspan)
	call getarg(12,c_T)
	call getarg(13,c_poly)	
	call getarg(14,c_npol)	
			
	read(c_cold,*) cold
	read(c_colv1,*) colv1     
	read(c_colv2,*) colv2     
	read(c_cole1,*) cole1     
	read(c_cole2,*) cole2     
	read(c_date_deb,*) date_deb
	read(c_date_fin,*) date_fin
	read(c_t0,*) t0	
	read(c_span,*)  span		        
	read(c_sspan,*) sspan	
	read(c_poly,*) poly
	read(c_npol,*) npol
	
	 
	fac = ISIGN(1,colv2) 
	colv2 = iabs(colv2) 
	     
c ########################       
c argument des harmoniques       
c ########################       
     	
	n = 0
	m = 1
	k = 1
	
	do while ( n.le. Len_trim(c_T) )
		n=n+1
		if ( c_T(n:n) .ne. '_' ) c_TT(k) = c_T(m:n)	
													
		if ( c_T(n:n) .eq. '_' ) then
			k = k + 1 			
			m = n+1			 			
		endif

		if ( n .eq. Len_trim(c_T) ) c_TT(k) = c_T(m:n)
		
	enddo 
	
        k1max = k 
	! write(6,*) 'Periodes selectionnees :' 
	do i =1,k1max
	read(c_TT(i),*) PERIOD(i) 
	! write (6,*) 'numero', i,':', PERIOD(i)	
	enddo
	
	if (PERIOD(1).eq.0) k1max=0
		                   	                                   
	! --------------------------       
	! lecture fichier excitation      
	! --------------------------             
       	maxc = max(cold,colv1)          
       	maxc = max(maxc,colv2)          
       	maxc = max(maxc,cole1)                      
       	maxc = max(maxc,cole2)
       
       	open(11,file=filein,status='old')   
        i = 0 
54     	read(11,*,end=65) (TAB(j),j=1,maxc)
       	i = i + 1 
	DAT1(i) = TAB(cold)              
	goto 54
65	continue
	close(11)
	
	if (date_deb.le.DAT1(1)) date_deb = DAT1(1)
	if (date_fin.ge.DAT1(i)) date_fin = DAT1(i)
	
	if (span.eq.0) span = date_fin - date_deb    
                                  
c --------------       
c Sliding window       
c --------------      
       
        open(20,file = '/home/www/html/eop-pc/prog_local/Temp/mca.res',status='unknown')        
        open(21,file = '/home/www/html/eop-pc/prog_local/Temp/mca_ipop.var',status ='unknown')       
        open(22,file = '/home/www/html/eop-pc/prog_local/Temp/mca_amph.var',status ='unknown')       

	
	iw  = 0       
4       iw  = iw + 1       
	
		date1 = date_deb + dble(iw-1) * sspan       
		date2 = date1 + span       
		datem = (date1 + date2) / 2.d0   ! date medianne   
      			          
       		
		open(11,file=filein,status='old')   

       		i = 0         
55     		read(11,*,end=66) (TAB(j),j=1,maxc)
       		date = TAB(cold)
       		xx   = TAB(colv1)             		
		if (cole1.ne.0) then
			sx   = TAB(cole1) 
		else
			sx = 1.d0
		endif

       		if (colv2.ne.0) then 
			yy  = TAB(colv2) * fac   			
       	              	if (cole2.ne.0) then
			sy  = TAB(cole2) 
				else
			sy = 1.d0
			endif
		endif
		
		if (date.lt.date1) goto 55        
		if (date.gt.date2) goto 66       
		i = i + 1
	 
        	DATG(i) = date     
        	 X(i)   = xx
		EX(i)   = sx	
       		if (colv2.ne.0) then 
			 Y(i) = yy     
            		EY(i) = sy
		endif
		
		
		goto 55      
      
 66		continue      
        	close(11)      
      
        	nmax = i       
   	   	 	    
      		do i = 1, nmax                             	
				
			dt = DATG(i) - t0 
		
			! --------------------------------        
			! observations matrix zp  et poids       
			! --------------------------------        
		  
			if (colv2.eq.0) then
				ZP(i) = X(i)  				    	        	  	       
				ZP(i) = 20 * cos(dpi/365.25* dt) + 10 * cos(dpi/433.0 * dt) 
				PP(i) = 1.d0 / ( EX(i)**2 +1d-8 )				
			else
				ZP(2*i-1) = X(i)  
				ZP(2*i  ) = Y(i)   
				!ZP(2*i-1) = 10 * dcos(dpi* dt) + 100 * dcos(dpi/1.2* dt) + 3 * dt +  5 * dcos(dpi/2.2* dt)
				!ZP(2*i  ) = 10 * dsin(dpi* dt) + 100 * dsin(dpi/1.2* dt)       	        	  	       
				PP(2*i-1) = 1.d0 / ( EX(i)**2 +1d-8 )
				PP(2*i  ) = 1.d0 / ( EY(i)**2 +1d-8 )			      
			endif
			 
			
			if (k1max.ne.0) then     		
				! argument a la date i
				do k = 1,k1max	
					A(k) = dpi / PERIOD(k) * dt   	
					A(k) = dmod(A(k),dpi) 			  
	  			enddo			          
				
				! ---------------------------------------------------        
				! observations relation matrix za : 
				!      x = \sum [m_c cos A(j) + m_s sin A(j)]             
				! ---------------------------------------------------            
      
     	        		if (colv2.eq.0) then
					do j =1,k1max       
	       					ZA(i,2*j-1) =  dcos(A(j))	
	        				ZA(i,2*j  ) =  dsin(A(j))
					enddo 
				else
				! ---------------------------------------------------
				! observations relation matrix za : 
				! 	x + i y = \sum  [ m_c  + i m_s ] exp (i a(j))	 
				! 	   chix = \sum [m_c cos A(j) - m_s sin A(j)]	   
				! 	   chiy = \sum [m_c sin A(j) + m_s cos A(j)]	   
				! ---------------------------------------------------            	
					do j =1,k1max       
	       					ZA(2*i-1,2*j-1) =  dcos(A(j))	
	        				ZA(2*i-1,2*j  ) = -dsin(A(j))
	       					ZA(2*i  ,2*j-1) =  dsin(A(j))	
	        				ZA(2*i  ,2*j  ) =  dcos(A(j))						  
					enddo            			
				endif
			endif 
				 
	        	
 			if (poly.eq.1) then 
				! estimation du polynome d'ordre n
				if (colv2.eq.0) then
					do j = 1, npol + 1      	  
	        				ZA(i,2*k1max +j) = dt**(j-1)      
	        			enddo
				else
					do j = 1, npol + 1      	  
	        				ZA(2*i-1,2*k1max + 2*j-1) = dt**(j-1)
						ZA(2*i-1,2*k1max + 2*j )  = 0.d0
						ZA(2*i  ,2*k1max + 2*j-1) = 0.d0
						ZA(2*i  ,2*k1max + 2*j  ) = dt**(j-1)	
	        			enddo
	     	        	endif
			endif	                 	       
        	enddo       
       
      		close(11)       
                     		
		if (colv2.eq.0) then
			n = nmax     ! dimension de la matrice colonne zp  
			if (poly.eq.0) l = 2 * k1max
			if (poly.eq.1) l = 2 * k1max + npol + 1      ! nombre de parametres a estimer     
        	else
			n = 2 * nmax ! dimension de la matrice colonne zp  
			if (poly.eq.0) l = 2 * k1max
			if (poly.eq.1) l = 2 * k1max + 2* (npol + 1) ! nombre de parametres a estimer     	
		endif  
			       
		! #######################        
		! least square treatement        
		! #######################          		          
				 
		call mcarre(n,l,za,zp,pp,xs,xr,xdx,cor,nes,rms,rms_r)        
		
		! #########################################   
		! Amplitudes et phase of the circular terms   
		! #########################################   
 	  
		if (c_t.ne.'0') then  
			do i = 1, k1max   
				call AMPHASE(XS(2*i-1),XDX(2*i-1),XS(2*i),XDX(2*i),amp,da,phase,dphase)   
        			AP(i)  = amp     
        			DAP(i) = da    
        			PH(i)  = phase   
        			DPH(i) = dphase 	 
				! procedure assurant la continuite de la phase d'une date a l'autre
				if (iw.ge.2) then 
				  if ( PH(i)-ph1.ge. 350.d0) PH(i) = PH(i) - 360.d0
				  if ( PH(i)-ph1.le.-350.d0) PH(i) = PH(i) + 360.d0	        			
				endif
				if (iw.eq.1) ph1 = PH(i)
			enddo   
   		endif 
		! ####################################################        
		! writing of the results in the file result polaam.res        
		! ####################################################        
        

      		write(20,1010) filein,date1,date2,n,l,t0,rms,rms_r,nes      
      		        
1010  		format ( 'INPUT FILE     --> ',a80,/,
     &                   'DATE DEP       --> ',f10.3,/,
     &		         'DATE END       --> ',f10.3,/,			  
     &		         'NB OBSERV      --> ',i10,/,
     &		         'NB PARAMETERS  --> ',i10,/,
     &			 'REF EPOCH      --> ',f10.3,/,
     &			 '                   ',/,
     &			 'RMS OBS        --> ',f10.3,/,
     &			 'RMS RESIDUAL   --> ',f10.3,/,
     &			 'NB REMOVED EQ. --> ',i10,/)
     
	       	    

1050  		format(
     &		'------------------------------------------------------',/,        
     &          '      x = m_c cos(2 pi/T * t) + m_s sin(2 pi/T * t)   ',/,
     &		'      T             m_c                    m_s        ',/,       
     &		'------------------------------------------------------')       
	 
	
1051  		format(
     &		'--------------------------------------------------------',/,        
     &          '           x + i y = (m_c  + i m_s ) exp (i 2 pi/T * t )',/,
     &		'       T             m_c                    m_s         ',/,       
     &		'--------------------------------------------------------')       
       	 
1052  		format(
     &		'--------------------------------------------------------',/,        
     &		'                    amplitude              phase        ',/,       
     &		'--------------------------------------------------------')       

1053  		format(
     &		'----------------------------',/,	 
     &          '           Polynomial term  ',/,
     &		'Power             x         ',/,	
     &		'----------------------------')       

1054  		format(
     &		'--------------------------------------------------------',/,        
     &          '                          Polynomial term               ',/,
     &		' Power               x                      y           ',/,       
     &		'--------------------------------------------------------')       

         	! Affichage des termes périodiques cos/sin
		if (colv2.eq.0) write(20,1050) 
		if (colv2.ne.0) write(20,1051)    
	
		if (k1max.ne.0) then
			do i = 1, k1max             	
	 			write(20,1061) PERIOD(i),XS(2*i-1),XDX(2*i-1),XS(2*i),XDX(2*i)     
         		enddo         
1061     		format( f9.2,4x ,f9.4,' +-',f8.4,3x,f9.4,' +-',f8.4,/)     
     		endif
 
         	! ------------------------------------------
		! Affichage des termes périodiques amp/phase
		if (colv2.ne.0) then
		write(20,1052)    
		if (k1max.ne.0) then
			do i = 1, k1max                	
	 			write(20,1062) PERIOD(i),AP(i),DAP(i),PH(i),DPH(i)     
         		enddo         
1062     		format( f9.2,4x ,f9.4,' +-',f8.4,3x,f9.1,' +-',f8.1,/)     
     		endif     
                endif  
                ! -------------------------------------------
		! Affichage des termes constants et seculaire  
                if (poly.eq.1) then
			if (colv2.eq.0)  write(20,1053)
			if (colv2.ne.0)  write(20,1054)
			do i = 1, npol+1 
 	 			if (colv2.eq.0) write(20,1063) i-1,xs(2*k1max+i),xdx(2*k1max+i)
				if (colv2.ge.1) write(20,1064) i-1,xs(2*k1max+2*i-1),xdx(2*k1max+2*i-1),xs(2*k1max+2*i),xdx(2*k1max+2*i)
 			enddo		
1063     		format('  ',i3,4x ,f9.4,' +-',f8.4,/)    
1064     		format('  ',i3,4x ,f9.4,' +-',f8.4,4x ,f9.4,' +-',f8.4,/)    
             	endif 
	 
		! -----------------------------------------
		! Sliding window :  in-phase / out-of-phase                
         	write(21,*) datem,(xs(2*i-1),xdx(2*i-1),xs(2*i),xdx(2*i),i=1,k1max)  	  
         	! Sliding window :  amplitude / phase 
		write(22,*) datem,(AP(i),DAP(i),PH(i),DPH(i),i=1,k1max)  	  
	          

		!#############       
		! correlations        
		!#############        
       
		write (20,2000)             
2000    	format(/,' Significant correlation coefficients > 0.1',/)       
       
    		do i=1,l       
			do j=1,l       
	  			if (i.eq.j) goto 3010       
	  			if ((cor(i,j).gt.0.1d0).and.(i.lt.j)) write(20,'(" cor(",i2,",",i2,")=",f5.2)') i,j,cor(i,j)             	    
3010    		enddo       
		enddo       
                            
		! ##########################        
		! residuals et modele estime           
		! ##########################                
       		open(18,file='/home/www/html/eop-pc/prog_local/Temp/mca.resi',status='unknown')        
        	do  i=1,nmax        
       			dt = DATG(i) - t0
			if (colv2.eq.0) then 
				X_mod = 0.d0
				Y_mod = 0.d0
				do k=1,k1max				
					A(k) = dpi / PERIOD(k) * dt 
					X_mod = X_mod + XS(2*k-1) * dcos(A(k)) +  XS(2*k) * dsin(A(k))  
				enddo
				
				do k=1,npol+1
				X_mod = X_mod + XS(2*k1max+k) * dt**(k-1) 
				enddo
			else
				X_mod = 0.d0
				Y_mod = 0.d0
				do k=1,k1max				
					A(k) = dpi / PERIOD(k) * dt 
					X_mod = X_mod + XS(2*k-1) * dcos(A(k)) -  XS(2*k) * dsin(A(k))   	
					Y_mod = Y_mod + XS(2*k-1) * dsin(A(k)) +  XS(2*k) * dcos(A(k))   
				enddo
				
				do k=1,npol+1
					X_mod = X_mod + XS(2*k1max+2*k-1) * dt**(k-1)
					Y_mod = Y_mod + XS(2*k1max+2*k)   * dt**(k-1) 					
				enddo				
			endif
			
			if (colv2.eq.0) write(18,2020) DATG(i),XR(i),X_mod  
       			if (colv2.ge.1) write(18,2020) DATG(i),XR(2*i-1),XR(2*i),X_mod,Y_mod  
					             
       		enddo        
       		close(18)        
2020   		format(5e20.8)
		

        	if (sspan.eq.0) goto 7 
		if (date2 + sspan.le.date_fin) goto 4     
  		if (date2 + sspan.gt.date_fin) goto 7    
  	          	        	  
 7      	continue  	
              
		close(21)       
		close(20)                           
    
	END        
       
       
c #################################################################       
        SUBROUTINE MCARRE(N,L,ZA,ZP,PP,XS,XR,XDX,COR,NES,rms,rms_r)       
c #################################################################       
           
        implicit double precision (a-h,o-z)             
	PARAMETER (n1=50000,k1=20,k2=2*k1+4)
       
        double precision ZA(n1,k2),ZP(n1),PP(n1),XR(n1),XS(k2),XDX(k2),COR(k2,k2),rms,rms_r   
        double precision ZA1(n,l),ZP1(n),ZS(n),ZB(l,l),ZC(l,l)      

        integer nopt,nes
	       
        nes  = 0                                                
        ! nopt = 1 ! option de supression d'equation d'observations
       	nopt = 0
	kt = 0 
 	! calcul de la variance des observations 
	xds = 0.d0        
	spp = 0.d0
        do i=1,n                	      
                spp = spp + PP(i)            ! somme des poids   
		xds = xds + ZP(i)**2 * PP(i) ! variance des observations		 	     
	enddo   
        rms = dsqrt( xds / spp )
	
	! ---------------------------------------------        
	! matrixs ZA1=ZA t(sqrt(ZP)) and ZP1=sqrt(ZP) L        
	! ---------------------------------------------         
 100    DO I=1,N        
        DO J=1,L        
        ZA1(I,J) = ZA(I,J)*DSQRT(PP(I)) 	  
        ENDDO        
        ZP1(I) = ZP(I)*DSQRT(PP(I))        	  
	ENDDO        
       
        do i=1,l
        ZS(I)=0.D0  
        enddo
       
        DO 24 I=1,l        
        DO 25 J=1,l        
        ZB(I,J)=0.D0        
 25     ZC(I,J)=0.D0        
 24     ZC(I,I)=1.D0        
        
        DO I=1,l       
        XS(I)=0.D0        
        ENDDO        
        
        ! ---------------------------- 
        ! normal equations : matrix ZB        
        ! ----------------------------         
        DO i=1,l                
        	DO j=1,l                
        		DO k=1,n        
        		ZB(i,j) = ZB(i,j) + ZA1(k,i) * ZA1(k,j)             
  			enddo
   		enddo       
        
        	DO k=1,n           
        	ZS(i) = ZS(i) + ZP1(k) * ZA1(k,i)        
		enddo
	        
   	enddo                
       
        ! ----------------        
        ! ZC: ZB INVERSION        
        ! ----------------        
        DO k=1,l                
        X=ZB(k,k)        
        	do J=1,L        
        	ZB(k,j) = ZB(k,j)/X                
        	ZC(k,j) = ZC(k,j)/X        
        	enddo	
        	DO 32 i=1,l                
        	if(i.eq.k) go to 32                
        	Y = ZB(i,k)                
        	DO j=1,l        
        	ZB(i,j) = ZB(i,j) - Y * ZB(k,j)		 
         	ZC(i,j) = ZC(i,j) - Y * ZC(k,j)	 
        	enddo		
32              CONTINUE                
   	enddo        
        
        ! -----------------------------        
        ! solutions XS and residuals XR        
        ! -----------------------------               
        DO i=1,l	
        DO j=1,l	
        XS(i)=XS(i)+ZS(j)*ZC(i,j)        
  	enddo
	enddo
	                
	! calcul de la variance des residus et du chi2
	xds = 0.D0        
	spp = 0.d0
        do i=1,n                	      
		XR(I) = ZP(i)                
        	do J=1,L        
        		XR(I) = XR(I) - XS(J) * ZA(I,J)        
 		enddo
                spp = spp + PP(i)            ! somme des poids   
		XDS = XDS + XR(i)**2 * PP(i) ! variance des residus		 
	enddo                                
        
	chi2 = xds / dble(n-l-kt)  ! chi2 
	
	rms_r  = dsqrt(xds / spp)  ! RMS RESIDUS        
     	
        ! ------------------------------------------       
        ! precision XDX and correlation coefficients        
        ! ------------------------------------------       
	        
        do i=1,l        
        	XDX(I)= dsqrt(chi2) * DSQRT(dabs(ZC(I,I)))   
	
        	do j=1,l        
        		cor(i,j) = dabs(zc(i,j) / dsqrt(DABS(zc(i,i)*zc(j,j))))        
		enddo
	enddo     
               
        ! ------------------        
        ! equations removing        
        ! ------------------         
	
	if (nopt.ne.0) then        
	    kt = 0              
      	  DO 52 i=1,n        
            XT  = 2.57d0 * rms_r               
	    XRM = DABS(XR(i))       
         
      	  IF(XRM.lT.XT) go to 52        
	    kt  = kt+1	   
	    nes = nes+1	     
         
	    do j=1,l	    
       	    ZA(i,j) = 0.d0	  
  	    enddo	
	    ZP(i) = 0.d0 
	    PP(i) = 0.d0       
 52 	continue  
       
	if (kt.ne.0) go to 100        
	         
	endif       
	       
	return              
    	end        
                   
C ########################################################################   
	SUBROUTINE AMPHASE(AR,DAR,AI,DAI,A,DA,PHASE,DPHASE)   
C COMPUTATION OF THE AMPLITUDE AND THE PHASE OF A COMPLEX NUMBER AR + I AI   
C       INPUT   :  AR,AI ....: COMPLEX NUMBER   
C		   DAR,DAI ..: ERRORS   
C	OUTPUT  :  A ........: AMPLITUDE   
C		   PHASE ....: PHASE   
C                  DA,DPHASE : ERRORS   
C #########################################################################   
	IMPLICIT REAL*8 (A-H,O-Z)    
   
	PI = 4.D0 * DATAN(1.D0)   
   
	A = DSQRT( AR**2 + AI**2)   
        
	DA = 1.D0/DSQRT(AR**2 + AI**2) * DABS(AR) * DAR +   
     &       1.D0/DSQRT(AR**2 + AI**2) * DABS(AI) * DAI    
   
           
	IF (AR.GT.0) PHASE = DATAN(AI/AR)   
           
        IF ((AR.EQ.0).AND.(AI.GT.0)) PHASE =  PI/2.D0    
        IF ((AR.EQ.0).AND.(AI.LT.0)) PHASE = -PI/2.D0    
   
  	IF (AR.LT.0) PHASE = DATAN(AI/AR) + PI   
   
	DPHASE = 1.D0/(1.D0 + (AI/AR)**2) * DAI /DABS(AR) +   
     &           1.D0/(1.D0 + (AI/AR)**2) * DAR * DABS(AI) /AR**2    
   
   
C CONVERSION IN DEGRES    
   	
	PHASE  = DMOD(PHASE  * 180.D0/PI,360.D0) 
        IF (PHASE.LT.0) PHASE = PHASE + 360.D0  
	DPHASE = DMOD(DPHASE * 180.D0/PI,360.D0)	   
	RETURN   
	END   
                                                                            
