program main
implicit none
!----------------
integer,parameter::nx=101,ny=51
real*8::p=0.3,q=0.3
real*8::dks=1.0,dat=1.0
real*8,dimension(nx,ny)::x=0,y=0,xp=0,yp=0
real*8,allocatable::xb0(:),yb0(:),xb1(:),yb1(:)
!---------------
integer::nb0,nb1
integer::i,ii,j,times,di,dj
!----------------
real*8::errmax=0.5,err
real*8::xks,xat,yks,yat,xksat,yksat,arfa,beta,gama,jj
real*8::a,ap,b,bp,c,cp,d,dp,e,ep
real*8::l,lp,m,mp,n,np
real*8::xe,xw,xn,xs,ye,yw,yn,ys
real*8::xx,yy
!---------------
real*8::thta,thta1,thta2,k1,k2,dis,dis1,disa,disb,disc
real*8::res
integer::iffind
!------------------------------------------------------------------------------------------------
!------------------------------------------------------------------------------------------------
!邊界坐標
open(1,file='boundary.txt') 
read(1,*) nb0 !左岸邊界
allocate(xb0(nb0),yb0(nb0))
read(1,*) ((xb0(i),yb0(i)),i=1,nb0)
read(1,*) nb1  !右岸岸邊界
allocate(xb1(nb1),yb1(nb1))
read(1,*) ((xb1(i),yb1(i)),i=1,nb1)
close(1)
!initial----------------------------------------------------------------------------------------
call cboundary(nx,nb0,xb0,yb0,x(:,1),y(:,1)) !(1,1) - (nx,1) 河床左岸邊界坐標
call cboundary(nx,nb1,xb1,yb1,x(:,ny),y(:,ny)) !(1,ny) - (nx,ny) 河床右岸邊界坐標
!------------------------------------
do i=1,nx,nx-1	!(1,1) - (1,ny),(nx,1) - (nx,ny) 河床上下邊界坐標
	do j=2,ny-1
		x(i,j)=x(i,1)+(x(i,ny)-x(i,1))*(j-1)/(ny-1)
		y(i,j)=y(i,1)+(y(i,ny)-y(i,1))*(j-1)/(ny-1)
	enddo
enddo
do j=2,ny-1 !(x,y) 河床內部節點邊界
	do i=2,nx-1
		x(i,j)=x(1,j)+(x(nx,j)-x(1,j))*(i-1)/(nx-1)
		y(i,j)=y(1,j)+(y(nx,j)-y(1,j))*(i-1)/(nx-1)
	enddo
enddo
do i=2,nx-1
	do j=2,ny-1
		x(i,j)=x(i,1)+(x(i,ny)-x(i,1))*(j-1)/(ny-1)
		y(i,j)=y(i,1)+(y(i,ny)-y(i,1))*(j-1)/(ny-1)
	enddo
enddo
!--------
open(1,file='initial.dat')
write(1,*) '"variable=x,y"'
write(1,*) 'zone  i=',ny,  'j=',nx
do i=1,nx
	do j=1,ny
		write(1,*) x(i,j),y(i,j)
	enddo
enddo
close(1)
!-------------------------------------------------------------------------------------------------
!-------------------------------------------------------------------------------------------------
times=0
100 continue
err=0
do i=2,nx-1
	do j=2,ny-1
		xks=(x(i+1,j)-x(i-1,j))/(2*dks)
		xat=(x(i,j+1)-x(i,j-1))/(2*dat)
		yks=(y(i+1,j)-y(i-1,j))/(2*dks)
		yat=(y(i,j+1)-y(i,j-1))/(2*dat)
		xksat=(x(i+1,j+1)+x(i-1,j-1)-x(i+1,j-1)-x(i-1,j+1))/(4*dks*dat)
		yksat=(y(i+1,j+1)+y(i-1,j-1)-y(i+1,j-1)-y(i-1,j+1))/(4*dks*dat)
		!----------
		arfa=xat*xat+yat*yat
		gama=xks*xks+yks*yks
		beta=xksat+yksat
		jj=xks*yat-xat*yks
!		考慮P、Q-----------------------------------------------------------------------
!		a=arfa*yks*yks/(xks*xks+yks*yks)
!		b=gama*yat*yat/(xat*xat+yat*yat)
!		c=-arfa*xks*yks/(xks*xks+yks*yks)
!		d=-gama*xat*yat/(xat*xat+yat*yat)
!		e=-2*beta*xksat
!		ap=c
!		bp=d
!		cp=arfa*xks*xks/(xks*xks+yks*yks)
!		dp=gama*xat*xat/(xat*xat+yat*yat)
!		ep=-2*beta*yksat
!		!--------------------------------
!		xe=x(i+1,j)
!		xw=x(i-1,j)
!		xn=x(i,j+1)
!		xs=x(i,j-1)
!		ye=y(i+1,j)
!		yw=y(i-1,j)
!		yn=y(i,j+1)
!		ys=y(i,j-1)
!		!--------------------------------
!		l=a*(xe+xw)/dks/dks+b*(xn+xs)/dat/dat+c*(ye+yw)/dks/dks+d*(yn+ys)/dat/dat+e
!		m=2*(a+b)
!		n=2*(c+d)
!		lp=ap*(xe+xw)/dks/dks+bp*(xn+xs)/dat/dat+cp*(ye+yw)/dks/dks+dp*(yn+ys)/dat/dat+ep
!		mp=2*(ap+bp)
!		np=2*(cp+dp)
!		!----------------------------
!		xp(i,j)=(l*np-lp*n)/(m*np-mp*n)
!		yp(i,j)=(l*mp-lp*m)/(mp*n-m*np)
!		if(abs(xp(i,j)-x(i,j))>err) err=abs(xp(i,j)-x(i,j))
!		if(abs(yp(i,j)-y(i,j))>err) err=abs(yp(i,j)-y(i,j))
!		x(i,j)=xp(i,j)
!		y(i,j)=yp(i,j)
!		不考慮P，Q----------------------------------------------------------------------
		p=0.
		q=0.
		!---------
		xp(i,j)=arfa*(x(i-1,j)+x(i+1,j))/dks/dks+gama*(x(i,j+1)+x(i,j-1))/dat/dat-2*beta*xksat+jj*(p*xks+q*xat)
		xp(i,j)=xp(i,j)/(2.0*(arfa/dks/dks+gama/dat/dat))
		if(abs(xp(i,j)-x(i,j))>err) err=abs(xp(i,j)-x(i,j))
		x(i,j)=xp(i,j)
		!---------
		yp(i,j)=arfa*(y(i-1,j)+y(i+1,j))/dks/dks+gama*(y(i,j+1)+y(i,j-1))/dat/dat-2*beta*yksat+jj*(p*yks+q*yat)
		yp(i,j)=yp(i,j)/(2.0*(arfa/dks/dks+gama/dat/dat))
		if(abs(yp(i,j)-y(i,j))>err) err=abs(yp(i,j)-y(i,j))
		y(i,j)=yp(i,j)
	enddo
enddo
!--------------------------------------------------------------------------------
goto 99
!邊界滑移------------------------------------------------------------------------
!x boundary----
!do j=1,ny,ny-1
!	if(j==1)  dj=1 !bottom
!	if(j==ny) dj=-1 !top
!	do i=2,nx-1
!		!print*,i,j
!		!尋找邊界控制點坐標
!		xx=0
!		yy=0
!		iffind=0
!		if(j==1)then
!			do ii=1,nb0-1
!				if(x(i,j)>=min(xb0(ii),xb0(ii+1)).and.x(i,j)<max(xb0(ii+1),xb0(ii)))then
!					xx=xb0(ii)
!					yy=yb0(ii)
!					iffind=1
!					cycle
!				endif
!			enddo
!		else
!			do ii=1,nb1-1
!				if(x(i,j)>=xb1(ii).and.x(i,j)<xb1(ii+1))then
!					xx=xb1(ii)
!					yy=yb1(ii)
!					iffind=1
!					cycle
!				endif
!			enddo
!		endif
!		if(iffind/=1) then
!			print*,'找不到邊界控制點坐標'
!			stop
!		endif
!		!-----------------------------------------------------------------
!		!										  x(i,j+dj)
!		!											/  *
!		!										/	 / |
!		!                             c     /       /  | 
!		!								/		 a /   | 
!		!						 	/			  /    | 
!		!						 /		 b		 /thta | 
!		!              ---------*---------------*--------------*----------
!		!					 (xx,yy)		  x(i,j)       (xx1,yy2)
!		!------------------------------------------------------------------
!		disa=sqrt((x(i,j+dj)-x(i,j))**2.0+(y(i,j+dj)-y(i,j))**2.0)		
!		disb=sqrt((xx-x(i,j))**2.0+(yy-y(i,j))**2.0)
!		disc=sqrt((x(i,j+dj)-xx)**2.0+(y(i,j+dj)-yy)**2.0)
!		if(disb/=0.and.disa/=0)then
!			thta=(disa*disa+disb*disb-disc*disc)/(2*disa*disb)
!			if(abs(thta)>1) thta=thta/abs(thta)
!			thta=acos(thta)
!			thta=3.1415926-thta
!			x(i,j)=x(i,j)+(x(i,j)-xx)*disa*cos(thta)/disb
!			y(i,j)=y(i,j)+(y(i,j)-yy)*disa*cos(thta)/disb
!		endif			
!!		disb=sqrt((x(i+di,j)-x(i,j))**2.0+(y(i+di,j)-y(i,j))**2.0)
!!		disc=sqrt((x(i,j+dj)-x(i+di,j))**2.0+(y(i,j+dj)-y(i+di,j))**2.0)
!!		thta=acos((disa*disa+disb*disb-disc*disc)/(2*disa*disb))
!!		x(i,j)=x(i,j)+(x(i+di,j)-x(i,j))*disa*cos(thta)/disb
!!		y(i,j)=y(i,j)+(y(i+di,j)-y(i,j))*disa*cos(thta)/disb
!	enddo
!enddo
!y boundary
!do i=1,nx,nx-1
!	if(i==1)  di=1
!	if(i==nx) di=-1	
!	do j=2,ny-1
!		disa=sqrt((x(i+di,j)-x(i,j))**2.0+(y(i+di,j)-y(i,j))**2.0)
!		if(j<ny)then
!			dj=1
!		elseif(j==ny)then
!			dj=-1
!		endif
!		disb=sqrt((x(i,j+dj)-x(i,j))**2.0+(y(i,j+dj)-y(i,j))**2.0)
!		disc=sqrt((x(i+di,j)-x(i,j+dj))**2.0+(y(i+di,j)-y(i,j+dj))**2.0)
!		thta=acos((disa*disa+disb*disb-disc*disc)/(2*disa*disb))
!		x(i,j)=x(i,j)+(x(i,j+dj)-x(i,j))*disa*cos(thta)/disb
!		y(i,j)=y(i,j)+(y(i,j+dj)-y(i,j))*disa*cos(thta)/disb
!	enddo
!enddo
!----------------------------------------------------
99 continue
!----------------------------------------------------
if(err>errmax)then
	times=times+1
	if(mod(times,500)==0)then
		print*,'times=',times,'   err=',err
		!if(times==10)then
		!	goto 101
		!endif
	endif
	goto 100
endif
!---------------------------------------------------------------------------
101 continue
!---------------------------------------------------------------------------
open(1,file='res.dat')
write(1,*) '"variable=x,y"'
write(1,*) 'zone  i=',ny,  'j=',nx
do i=1,nx
	do j=1,ny
		write(1,*) x(i,j),y(i,j)
	enddo
enddo
close(1)
!----------------------------------------------------------------------------
print*,'-------finished-------'
print*,'err=',err,' times=',times
!----------------------------------------------------------------------------
end




subroutine cboundary(nx,nb,xb,yb,x,y)
implicit none
!--------------------------------------------------
integer::nx,nb
real*8,dimension(nb)::xb,yb
real*8,dimension(nx)::x,y
!--------------------------------------------------
real*8::dis,dis1,disa,disb,disc
integer::i,ii
!-----------------------------------------------------------------------------------------------
!-----------------------------------------------------------------------------------------------
dis=0
do i=2,nb
	dis=dis+sqrt((xb(i)-xb(i-1))**2.0+(yb(i)-yb(i-1))**2.0)
enddo
!---------------------------------------------------
x(1)=xb(1)
y(1)=yb(1)
x(nx)=xb(nb)
y(nx)=yb(nb)
dis1=dis/(nx-1)
disa=0
!-------------------
do i=2,nx-1
	disa=(i-1)*dis1
	disb=sqrt((xb(1)-xb(2))**2.0+(yb(1)-yb(2))**2.0)
	if(disa<=disb)then
		x(i)=xb(1)+(xb(2)-xb(1))*disa/disb
		y(i)=yb(1)+(yb(2)-yb(1))*disa/disb
	else
		disb=0
		do ii=2,nb-1
			disb=disb+sqrt((xb(ii)-xb(ii-1))**2.0+(yb(ii)-yb(ii-1))**2.0)
			disc=disb+sqrt((xb(ii+1)-xb(ii))**2.0+(yb(ii+1)-yb(ii))**2.0)
			if(disa>disb.and.disa<=disc)then
				x(i)=xb(ii)+(xb(ii+1)-xb(ii))*(disa-disb)/(disc-disb)
				y(i)=yb(ii)+(yb(ii+1)-yb(ii))*(disa-disb)/(disc-disb)
				exit
			endif			
		enddo
		if(disa>disc)then
			print*,'cannot find position'
			print*,i,disa,disb,disc
			stop
		endif
	endif
enddo
!----------------------------------------------------------------------------------------------
return
end subroutine cboundary