program hw
implicit none
real,allocatable::x(:,:),y(:,:)
real::d,teta1,teta2,r2,r1,l,teta
INTEGER::np,nr,i,j
R1=6.0
R2=20.0
teta1=135.0
teta2=45.0
nr=40
np=50
ALLOCATE(x(np,nr),y(np,nr))

!Botttom line Grid #1
call grider(-R1*sind(teta1-90),R1*sind(90-teta2),np,X(1:np,1))
Y(1:np,1)=R1*cosd(teta1-90)

!Top line Grid #2
call grider(teta1,teta2,np,X(1:np,nr))
Y(1:np,nr)=R2*sind(X(1:np,nr))
X(1:np,nr)=R2*cosd(X(1:np,nr))

!middle lines Grid #3
do i=1,np
call grider(x(i,1),x(i,nr),nr,x(i,1:nr))
call grider(y(i,1),y(i,nr),nr,y(i,1:nr))
end do

!plot #4
open(unit=100,file='mesh.plt')
write(100,*)"variables = x ,y"
write(100,*)"ZONE I=",nr,"J=",np
do i=1,np
        do j=1,nr
        write(100,*)X(i,j),Y(i,j)   
        end do
end do
close(100)
end program

SUBROUTINE grider(L1,L2,n,y)
real::etha,beta,alpha,y(n),L1,L2
INTEGER::n,i
alpha=0.5
beta=1.1
y=1
Do i=1,n
        etha = real(i-1)/real(n-1)
        y(i) = L1+real(L2-L1) * ((2*alpha + beta) * (((beta+1)/(beta-1))**((etha-alpha)/(1-alpha))) + 2*alpha - beta) /&
         ((2*alpha+1) * (1+((beta+1)/(beta-1))**((etha-alpha)/(1-alpha))))
end do
end SUBROUTINE