c+++
c
c	PROGRAM		WAVINESS_GEN
c
c	PURPOSE		This program generates a random slope error surface
c
c	ALGORITHM       See M. Sanchez del Rio and A. Marcelli,
c			NIM-A319 (1992) 170-177
c
c	COMPILATION NOTE: f77 +E1 in hp systeme (MSR 95/09/25)
c
c	MODIFICATION HISTORY
c			June 1991	first version
c			September 1991	update 
c			October 1991	moved to Sun
c			September 1992	distribution version
c---
	program		waviness_gen
	implicit real*4 (a-h,o-z)
	character*80	file
	real*4	x(101),y(101),coe(100),shi(100),ranc(100),shin(100)
	real*4	zz(201,201),xx(201),yy(201)
	data	pi	/3.1416/
	data	twopi /6.2832/
c
c  reads input
c
	write (6,*) 'Name of output file for shadow ? '
	read (*,3333) file
3333	format (a)
	open (21,file=file,status='unknown')
c
	write (6,*) ' input number of points along the X(transversal) and'
	write (6,*) ' Y(along) directions of the mirror ? '
	read(*,*) npointx,npointy
	write (6,*) 
     $'input mirror width and mirror length [in user units] ? '
	read(*,*) width,xlength
	write (6,*) ' input # of harmonics to be considered and an
     $ estimation of the wanted slope error [arc sec] ? '
	read(*,*) n,slp
	yfact = 1.
	slp = slp*2.777e-4*pi/180.   ! pass to rads
c
	xfact = xlength/pi
	i_plot = 0
	xntot = 0
	write (6,*) ' Input the seed for the random number generator : '
	read(*,*) iseed
	write (6,*) 'Now for each harmonic input three values:        '
	write (6,*) ' i)   C:  weighting coefficient  0<=C<=1         '
	write (6,*) ' ii)  y: an initial random phase shift           '
	write (6,*) ' iii) g: a fraction of the total lenght interval '
	do 1011 i=1,n+1	
	 write (6,*) ' for harmonic #  ',i-1,' input  C,y,g ? '
	 read(*,*) a1,a2,a3
	 xntot = xntot + a1
	 amp = slp*xlength/(0.5+float(i-1))/sqrt(2.)/pi
	 coe(i) = a1*amp
	 shin(i) = a2
	 ranc(i) = a3 
	 write (6,*) coe(i),shin(i),ranc(i)
1011	continue

c	write (6,*) 'xntot = ',xntot
	do 10111 i=1,n+1
	coe(i) = coe(i)/sqrt(xntot)       ! for a single harmonic
10111	continue
c
c	begin calculations
c
	stp = pi/(npointy-1)
	stpx = width/(npointx-1)
	stpy = xlength/(npointy-1)
	k=0
10 	continue
	k=k+1
	xx(k)=-(width/2)+stpx*(k-1)
	do 1022 j=1,n+1
c	tmp = rand(iseed)
	tmp = ran(iseed)
	shi(j)=shin(j)*pi+tmp*pi*ranc(j)
1022	continue

	do 1033 i=1,npointy
	y(i)=0
1033	continue

	do 1044 j=1,n+1
	do 1055 i=1,npointy
	 x(i)= -pi/2 + (i-1)*stp
	 nn=float(j-1)
c	 ynew= coe(j)*4/pi*( (-1)**nn*cos((2*nn+1)*(x(i)+shi(j)))/(2*nn+1))
	 ynew= coe(j)*( (-1)**nn*cos((2*nn+1)*(x(i)+shi(j))))
	 y(i)=ynew+y(i)
1055	continue
1044	continue
	do 1066 i=1,npointy
	 x(i)=x(i)*xfact
	 yy(i) = x(i)
	 y(i)=y(i)*yfact
	 zz(k,i)=y(i)
1066	continue
	if (k.lt.npointx) go to 10
c
c	writes file for shadow
c
	write (21,*) npointx,npointy
	write (21,*) (yy(i), i=1,npointy)
	do 1099 i=1,npointx
	write (21,*) xx(i),(zz(i,k),k=1,npointy)
1099	continue
	close(21)
c
c  calculate the slope error rms of the whole surface
c
	call surrms(npointy,npointx,stpy,stpx,zz,stdy,stdx)
	write (6,*) ' '
	write (6,*) '	Slope error RMS along Y is: ',stdy,' arc sec'
	write (6,*) '	Slope error RMS along X is: ',stdx,' arc sec'
	write (6,*) ' '
	write (6,*) '   Output file is : ',file
	write (6,*) '   Don`t forget to run PRESURFACE before SHADOW'
	write (6,*) ' '
c
	end
C+++
C
C
C	subroutine surrms: calculates the slope error RMS
C	of a surface in parallel and perpendicular direction
C
C
C
C---
	subroutine surrms(imax,jmax,d,dx,z,stdy,stdx)
	real*4		diff(201,201)
	real*4		z(201,201)
*calculates RMS along y
	difftot=0
	do j=1,jmax
	 do i=1,imax-1
	  diff_d=z(j,i+1)-z(j,i)
	  diff(j,i)=(57.29578/2.777e-4)*atan(diff_d/d)
	  difftot=difftot+diff(j,i)
	 end do
	end do	
	do j=jmax+1,201
	 do i=imax,201
	  diff(j,i)=0.0
	 end do
	end do	
*
	ndiff=(jmax*(imax-1))
	diffmean=difftot/ndiff
	diffstd = 0.0
	do j=1,jmax
	 do i=1,imax-1
	  diffstd = diffstd+(diff(j,i)-diffmean)**2
	 end do
	end do	
	stdy = sqrt(diffstd/(ndiff-2))		
*calculates RMS along x
	difftot=0
	do i=1,imax
	 do j=1,jmax-1
	  diff_d=z(j,i)-z(j+1,i)
	  diff(j,i)=(57.29578/2.777e-4)*atan(diff_d/dx)
	  difftot=difftot+diff(j,i)
	 end do
	end do	
	do j=jmax,201
	 do i=imax+1,201
	  diff(j,i)=0.0
	 end do
	end do	
*
	ndiff=((jmax-1)*imax)
	diffmean=difftot/ndiff
	diffstd = 0.0
	do j=1,jmax-1
	 do i=1,imax
	  diffstd = diffstd+(diff(j,i)-diffmean)**2
	 end do
	end do	
	write (*,*) '>>>>>>> diffstd: ',diffstd
	stdx = sqrt(diffstd/(ndiff-2))		
	write (*,*) '>>>>>>> stdx: ',stdx
c
	return 
	end