!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
! |
! ||===\\ |
! || \\ |
! || || //==\\ || || //==|| ||/==\\ |
! || || || || || || || || || || |
! || // || || || || || || || |
! ||===// \\==// \\==\\ \\==\\ || |
! |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
! |
! CREATE_SURF May. 2008 |
! |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
subroutine create_surf(surface,is,debug)
!------------------------------------------------------------------------------|
!(((((((((((((((( Purpose of the routine ))))))))))))))))))))))))))))))))))))))
!------------------------------------------------------------------------------|
! for a given surfaces passed as argument, this routine computes
! the 8 fields x,y,z,xn,yn,zn,r,s for each point of the surface
! according to its type and levelt, and the triangulation
!------------------------------------------------------------------------------|
!(((((((((((((((( declaration of the subroutine arguments ))))))))))))))))))))
!------------------------------------------------------------------------------|
use constants
use definitions
!use mpi
use randomodule
implicit none
include 'mpif.h'
type (sheet) :: surface
integer is
integer debug
!------------------------------------------------------------------------------|
!(((((((((((((((( declaration of the subroutine internal variables )))))))))))))
!------------------------------------------------------------------------------|
integer :: levelt,surface_type,nproc,iproc,ierr,i,j,seed,indix
integer :: err,ntmax,nhmax,npmax,nmax,nmode,nvmax,nnpnmax,nh,nohalt_hull,loc
integer :: dmode,nt,nhull,icircles,irect
integer,dimension(:), allocatable :: hulltriangles,vis_tlist,vis_elist
integer,dimension(:), allocatable :: add_tlist
integer,dimension(:,:),allocatable :: vertices,neighbour
logical*1,dimension(:),allocatable :: lt_work,ln_work
double precision :: sp01,sp02,sp03,sp04,sp05,sp06,sp07,sp08,sp09,sp10,sp11,sp12
double precision :: sp13,sp14
double precision :: delta,eps,epsil,delta_angle,delta_radius,angle
double precision :: xx,yy,xini,xend,yini,yend,dx,dy,dist
double precision,dimension(:), allocatable :: field
double precision,dimension(:,:),allocatable :: points,centres
logical clockwise
character ch
character*72 :: shift
integer nnn(1),nnlist(1),ntrilist(1)
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
call mpi_comm_size (mpi_comm_world,nproc,ierr)
call mpi_comm_rank (mpi_comm_world,iproc,ierr)
shift=' '
if (surface%surface_type.lt.0) then
!---------------------------------
! reading surface from ascii file
!---------------------------------
if (iproc.eq.0) then
write (ch,'(i1)') -surface%surface_type
write(*,*) shift//'surface-'//ch//'.txt'
open(71,file='surface-'//ch//'.txt',status='old')
read(71,*) surface%nsurface,surface%nt
write(*,*) shift//'-> nb of points :',surface%nsurface
write(*,*) shift//'-> nb of triangles:',surface%nt
end if
call mpi_bcast (surface%nsurface,1,mpi_integer,0,mpi_comm_world,ierr)
call mpi_bcast (surface%nt,1,mpi_integer,0,mpi_comm_world,ierr)
deallocate (surface%x,surface%y,surface%z)
deallocate (surface%xn,surface%yn,surface%zn)
deallocate (surface%r,surface%s)
deallocate (surface%icon)
allocate (surface%x(surface%nsurface),surface%y(surface%nsurface),surface%z(surface%nsurface))
allocate (surface%xn(surface%nsurface),surface%yn(surface%nsurface),surface%zn(surface%nsurface))
allocate (surface%r(surface%nsurface),surface%s(surface%nsurface))
allocate (surface%icon(3,surface%nt))
if (iproc.eq.0) then
write(*,*) surface%nsurface,surface%nt
do i=1,surface%nsurface
read (71,*) surface%x(i),surface%y(i),surface%z(i),surface%xn(i),surface%yn(i),surface%zn(i)
surface%r(i)=surface%x(i)
surface%s(i)=surface%y(i)
if (surface%x(i)<0 .or. surface%x(i)>1 .or. surface%y(i)<0 .or. surface%y(i)>1 .or. surface%z(i)<0 .or. surface%z(i)>1 &
.and. iproc.eq.0) &
print *,surface%x(i),surface%y(i),surface%z(i)
enddo
do i=1,surface%nt
read (71,*) (surface%icon(j,i),j=1,3)
enddo
close (71)
endif
call mpi_bcast(surface%x,surface%nsurface,mpi_double_precision,0,mpi_comm_world,ierr)
call mpi_bcast(surface%y,surface%nsurface,mpi_double_precision,0,mpi_comm_world,ierr)
call mpi_bcast(surface%z,surface%nsurface,mpi_double_precision,0,mpi_comm_world,ierr)
call mpi_bcast(surface%xn,surface%nsurface,mpi_double_precision,0,mpi_comm_world,ierr)
call mpi_bcast(surface%yn,surface%nsurface,mpi_double_precision,0,mpi_comm_world,ierr)
call mpi_bcast(surface%zn,surface%nsurface,mpi_double_precision,0,mpi_comm_world,ierr)
call mpi_bcast(surface%r,surface%nsurface,mpi_double_precision,0,mpi_comm_world,ierr)
call mpi_bcast(surface%s,surface%nsurface,mpi_double_precision,0,mpi_comm_world,ierr)
call mpi_bcast(surface%icon,surface%nt*3,mpi_integer,0,mpi_comm_world,ierr)
else
!----------------
! create surface
!----------------
seed=1234567
sp01=surface%sp01
sp02=surface%sp02
sp03=surface%sp03
sp04=surface%sp04
sp05=surface%sp05
sp06=surface%sp06
sp07=surface%sp07
sp08=surface%sp08
sp09=surface%sp09
sp10=surface%sp10
sp11=surface%sp11
sp12=surface%sp12
sp13=surface%sp13
sp14=surface%sp14
surface_type=surface%surface_type
levelt=surface%levelt
! delta=2.d0**(-levelt) ! spacing between points on hull
delta=1.d0 / (2.d0**levelt +1) ! opla
epsil=(surface%stretch-1)/2.d0
if (iproc.eq.0) write(8,*) 'create_surfaces: delta=',delta
if (iproc.eq.0) write(8,*) 'create_surfaces: epsilon=',epsil
!-----building the surface-------------
! type 1 : flat surface
! type 2 : rectangular emboss
! type 3 : convex spherical emboss
! type 4 : concave spherical emboss
! type 5 : double rectangular emboss
! type 6 : sinus
! type 7 : noisy surface
! type 8 : double sinus
! type 9 : cosinus sinus
! type 10 : slope
! type 11 : 2D hill
! type 12 : rectangular emboss with specified slope
! type 13 : slope with 2 kinks
select case(surface_type)
case (1,2,5,6,7,8,9,10,11,12,13)
! if (surface%rand) then
indix=1
! convex hull
do i=1,2**levelt +1
surface%x(indix)=dble(i-1)*delta
surface%y(indix)=0.d0
indix=indix+1
end do
do i=1,2**levelt +1
surface%x(indix)=1.d0
surface%y(indix)=dble(i-1)*delta
indix=indix+1
end do
do i=1,2**levelt +1
surface%x(indix)=1.d0-dble(i-1)*delta
surface%y(indix)=1.d0
indix=indix+1
end do
do i=1,2**levelt +1
surface%x(indix)=0.d0
surface%y(indix)=1.d0-dble(i-1)*delta
indix=indix+1
end do
!interior nodes
if(surface%rand) then
call random_number (surface%x(indix:surface%nsurface))
call random_number (surface%y(indix:surface%nsurface))
else
surface%x(indix:surface%nsurface)=0.5d0
surface%y(indix:surface%nsurface)=0.5d0
end if
do j=2,2**levelt +1
do i=2,2**levelt +1
surface%x(indix)=(((surface%x(indix)-0.5d0)*2.d0)*epsil+dble(i-1))*delta
surface%y(indix)=(((surface%y(indix)-0.5d0)*2.d0)*epsil+dble(j-1))*delta
indix=indix+1
enddo
enddo
! else
! ! regular grid of points
! indix=1
! do j=1,2**levelt+1
! do i=1,2**levelt+1
! surface%x(indix)=float(i-1)*delta
! surface%y(indix)=float(j-1)*delta
! indix=indix+1
! end do
! end do
! end if
case (22)
if (surface%rand) then
indix=1
! convex hull
do i=1,2**levelt
surface%x(indix)=(i-1)*delta
surface%y(indix)=0.d0
indix=indix+1
end do
do i=1,2**levelt
surface%x(indix)=1.d0
surface%y(indix)=(i-1)*delta
indix=indix+1
end do
do i=1,2**levelt
surface%x(indix)=1.d0-(i-1)*delta
surface%y(indix)=1.d0
indix=indix+1
end do
do i=1,2**levelt
surface%x(indix)=0.d0
surface%y(indix)=1.d0-(i-1)*delta
indix=indix+1
end do
if (iproc.eq.0) write(8,*) 'nhull=',indix-1,'for surface ',is
! placing rows of points outside, on and outside the rectangle
eps=0.001
do irect=-5,3,2
xini=sp02+irect*eps
xend=sp03-irect*eps
yini=sp04+irect*eps
yend=sp05-irect*eps
dx=(xend-xini)/2**levelt
dy=(yend-yini)/2**levelt
do i=1,2**levelt
xx=xini+(i-1)*dx
yy=yini
if (xx.ge.0 .and. xx.le.1 .and. yy.ge.0 .and. yy.le.1) then
surface%x(indix)=xx
surface%y(indix)=yy
indix=indix+1
end if
end do
do i=1,2**levelt
xx=xend
yy=yini+(i-1)*dy
if (xx.ge.0 .and. xx.le.1 .and. yy.ge.0 .and. yy.le.1) then
surface%x(indix)=xx
surface%y(indix)=yy
indix=indix+1
end if
end do
do i=1,2**levelt
xx=xend-(i-1)*dx
yy=yend
if (xx.ge.0 .and. xx.le.1 .and. yy.ge.0 .and. yy.le.1) then
surface%x(indix)=xx
surface%y(indix)=yy
indix=indix+1
end if
end do
do i=1,2**levelt
xx=xini
yy=yend-(i-1)*dy
if (xx.ge.0 .and. xx.le.1 .and. yy.ge.0 .and. yy.le.1) then
surface%x(indix)=xx
surface%y(indix)=yy
indix=indix+1
end if
end do
end do
! throwing remaining points at random
do while (indix.le.surface%nsurface)
xx=ran(seed)
yy=ran(seed)
if ((xx<sp02-10*eps .or. xx>sp03+10*eps .or. yy<sp04-10*eps .or. yy>sp05+10*eps) .or. &
(xx>sp02+10*eps .and. xx<sp03-10*eps .and. yy>sp04+10*eps .and. yy<sp05-10*eps)) then
surface%x(indix)=xx
surface%y(indix)=yy
indix=indix+1
end if
end do
else
! regular grid of points
indix=1
do j=1,2**levelt+1
do i=1,2**levelt+1
surface%x(indix)=float(i-1)*delta
surface%y(indix)=float(j-1)*delta
indix=indix+1
end do
end do
end if
case (3,4)
if (surface%rand) then
! convex hull
indix=1
do i=1,2**levelt
surface%x(indix)=(i-1)*delta
surface%y(indix)=0.d0
indix=indix+1
end do
do i=1,2**levelt
surface%x(indix)=1.d0
surface%y(indix)=(i-1)*delta
indix=indix+1
end do
do i=1,2**levelt
surface%x(indix)=1.d0-(i-1)*delta
surface%y(indix)=1.d0
indix=indix+1
end do
do i=1,2**levelt
surface%x(indix)=0.d0
surface%y(indix)=1.d0-(i-1)*delta
indix=indix+1
end do
nhull=indix-1
if (iproc.eq.0) write(8,*) 'nhull=',nhull,'on surface ',is
! placing rows of nhull points outside, in, and inside the circle
delta_angle=2.d0*pi/nhull
delta_radius=0.0025d0*sp04
do icircles=-4,2
do i=0,nhull-1
angle=i*delta_angle
surface%x(indix)=sp02+(sp04+icircles*delta_radius)*cos(angle)
surface%y(indix)=sp03+(sp04+icircles*delta_radius)*sin(angle)
indix=indix+1
end do
end do
! throwing the rest of the points at random provided they are not
! too close to the concentric circles so that it does not disturb
! the regular triangulation
do while (indix.le.surface%nsurface)
xx=ran(seed)
yy=ran(seed)
dist=sqrt((xx-sp02)**2+(yy-sp03)**2)
if (dist>(sp04+10*delta_radius) .or. dist<(sp04-10*delta_radius)) then
surface%x(indix)=xx
surface%y(indix)=yy
indix=indix+1
end if
end do
else
! regular grid of points
indix=1
do j=1,2**levelt+1
do i=1,2**levelt+1
surface%x(indix)=float(i-1)*delta
surface%y(indix)=float(j-1)*delta
indix=indix+1
end do
end do
end if
case default
call stop_run ('surface type not defined$')
end select
!-------------------------------------------|
!-----building icon array-------------------|
!-------------------------------------------|
! we use the NN library to compute the triangulation
! of points for the surface
ntmax=surface%nsurface*3
nhmax=surface%nsurface
npmax=surface%nsurface
nnpnmax=100
nmax=3*ntmax+npmax
nvmax=ntmax
allocate (points(2,surface%nsurface),stat=err) ; if (err.ne.0) call stop_run ('Error alloc points in create_surface$')
allocate (field(surface%nsurface),stat=err) ; if (err.ne.0) call stop_run ('Error alloc field in create_surface$')
allocate (vertices(3,ntmax),stat=err) ; if (err.ne.0) call stop_run ('Error alloc vertices in create_surface$')
allocate (centres(3,ntmax),stat=err) ; if (err.ne.0) call stop_run ('Error alloc centress in create_surface$')
allocate (neighbour(3,ntmax),stat=err) ; if (err.ne.0) call stop_run ('Error alloc neighbour in create_surface$')
allocate (hulltriangles(nhmax),stat=err) ; if (err.ne.0) call stop_run ('Error alloc hulltriangles in create_surface$')
allocate (vis_tlist(nvmax),stat=err) ; if (err.ne.0) call stop_run ('Error alloc vis_tlist in create_surface$')
allocate (vis_elist(nvmax),stat=err) ; if (err.ne.0) call stop_run ('Error alloc vis_elist in create_surface$')
allocate (add_tlist(nvmax),stat=err) ; if (err.ne.0) call stop_run ('Error alloc add_tlist in icreate_surface$')
allocate (lt_work(ntmax),stat=err) ; if (err.ne.0) call stop_run ('Error alloc lt_work in create_surface$')
allocate (ln_work(surface%nsurface),stat=err) ; if (err.ne.0) call stop_run ('Error alloc ln_work in create_surface$')
points(1,:)=surface%x
points(2,:)=surface%y
dmode=-2
nmode=0
loc=1
nohalt_hull=0
clockwise=.true.
eps=1.d-8
call nn2d_setup (surface%nsurface,ntmax,nhmax,npmax,nnpnmax,nmax, &
points,dmode,nmode,clockwise,field,nt,vertices, &
centres,neighbour,nh,hulltriangles,nohalt_hull, &
loc,nnn,nnlist,ntrilist, &
eps,nvmax,vis_tlist,vis_elist,add_tlist, &
lt_work,ln_work)
if (iproc.eq.0) write (8,*) 'surface ',is,' has ',nt,' triangles'
surface%icon(1:3,1:surface%nt)=vertices(1:3,1:nt)
deallocate (points)
deallocate (field)
deallocate (vertices)
deallocate (centres)
deallocate (neighbour)
deallocate (hulltriangles)
deallocate (vis_tlist)
deallocate (vis_elist)
deallocate (add_tlist)
deallocate (lt_work)
deallocate (ln_work)
!----------------------------------------------|
!----computing z coordinate--------------------|
!----------------------------------------------|
call zpoints (surface%nsurface,surface%x,surface%y,surface%z,surface_type, &
sp01,sp02,sp03,sp04,sp05,sp06,sp07,sp08,sp09,sp10,sp11,sp12, &
sp13,sp14)
!----------------------------------------------|
!----computing the normals---------------------|
!----------------------------------------------|
call compute_normals (surface%nsurface,surface%x,surface%y,surface%z, &
surface%nt,surface%icon,surface%xn,surface%yn,surface%zn)
if (debug>=2) call output_surf (surface,is,'init',-1,-1)
end if
end subroutine create_surf
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
subroutine zpoints (ns,x,y,z,surface_type,sp01,sp02,sp03,sp04,sp05,sp06,sp07, &
sp08,sp09,sp10,sp11,sp12,sp13,sp14)
!------------------------------------------------------------------------------|
!(((((((((((((((( Purpose of the routine ))))))))))))))))))))))))))))))))))))))
!------------------------------------------------------------------------------|
! knowing the two-dimensional spatial distribution of points on the plane (x,y)
! this routine computes the z-coordinate accordingly to the surface type.
! ns is the number of points
! x,y,z are the coordinates of the points
! surface_type is the type of surface under consideration
! sp01..sp14 are the surface parameters (not all are necessarily used)
!------------------------------------------------------------------------------|
!(((((((((((((((( declaration of the subroutine arguments ))))))))))))))))))))
!------------------------------------------------------------------------------|
use constants
implicit none
integer ns
double precision x(ns),y(ns),z(ns)
integer surface_type
double precision sp01,sp02,sp03,sp04,sp05,sp06,sp07,sp08,sp09,sp10,sp11,sp12
double precision sp13,sp14
!------------------------------------------------------------------------------|
!(((((((((((((((( declaration of the subroutine internal variables )))))))))))))
!------------------------------------------------------------------------------|
integer i
double precision dist,ran,m,x1a,x1b,x2a,x2b,y1a,y1b,y2a,y2b,zx
double precision psi,theta1,theta2,theta1s,theta2s,psis,y0s,y1s,y2s,x0,y0w
double precision y0a,y0b,y0,y1,y2,a1,b1,bb1,c1,m1,a2,b2,bb2,c2,m2,a3,b3,bb3,c3
double precision m3,a4,b4,bb4,c4,m4,eps,psi2,psi2s,psi3,psi3s,dz1,dz2,dz3,dz4
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
eps=1.d-10
select case(surface_type)
case (1)
! a flat surface,
! sp01 is the z level
z=sp01
case (2)
! rectangular emboss,
! sp01 is the z level
! sp02 and 03 are x1,x2
! sp04 and 05 are y1,y2
! sp06 is the thickness
do i=1,ns
if ((x(i)-sp02)*(x(i)-sp03).le.0.d0 .and. (y(i)-sp04)*(y(i)-sp05).le.0.d0) then
z(i)=sp01-sp06
else
z(i)=sp01
end if
end do
case (3)
! convex spherical emboss,
! sp01 is the z level
! sp02 and 03 are x0,y0
! sp04 is the radius
do i=1,ns
dist=sqrt((x(i)-sp02)**2+(y(i)-sp03)**2)
if (dist.le.sp04) then
z(i)=sp01+sp04*sin(acos(dist/sp04))
else
z(i)=sp01
endif
end do
case (4)
! convex spherical emboss,
! sp01 is the z level
! sp02 and 03 are x0,y0
! sp04 is the radius
do i=1,ns
dist=sqrt((x(i)-sp02)**2+(y(i)-sp03)**2)
if (dist.le.sp04) then
z(i)=sp01-sp04*sin(acos(dist/sp04))
else
z(i)=sp01
endif
end do
case (5)
! double rectangular emboss,
! sp01 is the z level
! sp02 and 03 are x1,x2
! sp04 and 05 are x3,x4
! sp06 and 07 are y1,y2
! sp08 and 09 are y3,y4
! sp10 is the thickness
z=sp01
do i=1,ns
if ((x(i)-sp02)*(x(i)-sp03).le.0.d0 .and. (y(i)-sp06)*(y(i)-sp07).le.0.d0) then
z(i)=z(i)-sp10
endif
if ((x(i)-sp04)*(x(i)-sp05).le.0.d0 .and. (y(i)-sp08)*(y(i)-sp09).le.0.d0) then
z(i)=z(i)-sp10
endif
end do
case (6)
! a sinus,
! sp01 is the z level
! sp02 is the wavelength
! sp03 is the amplitude
! sp04 is the phase
do i=1,ns
z(i)=sp01+sp03*sin(2.d0*pi*(x(i)/sp02+sp04))
end do
case (7)
! a noisy surface,
! sp01 is the z level
! sp02 is the noise amplitude
do i=1,ns
call random_number (ran)
z(i)=sp01+ran*sp02
end do
case (8)
! a double sinus.
! sp01 is the z level
! sp02 is the x-wavelength
! sp03 is the x-amplitude
! sp04 is the y-wavelength
! sp05 is the y-amplitude
do i=1,ns
z(i)=sp01+sp03*sin(x(i)*2.*pi/sp02)+sp05*sin(y(i)*2.*pi/sp04)
end do
case (9)
! a cosinus,
! sp01 is the z level
! sp02 is the wavelength
! sp03 is the amplitude
! sp04 is the phase
do i=1,ns
z(i)=sp01+sp03*cos(2.d0*pi*(x(i)/sp02+sp04))
end do
case (10)
! a 2D embankment
! sp01 is z0
! sp02 is y0
! sp03 is psi
! sp04 is the thickness
do i=1,ns
if (y(i)<sp02) then
z(i)=sp01
else
z(i)=sp01+min((y(i)-sp02)*tan(sp03*pi/180.d0),sp04)
end if
end do
case (11)
! a 2D hill
! sp01 is z level
! sp02 height
! sp03 is y0
! sp04 is width
do i=1,ns
if (y(i)>sp03-sp04/2.d0 .and. y(i)<sp03+sp04/2.d0) then
z(i)=sp01+sp02*(sin(pi*((y(i)-sp03)/sp04+0.5d0)))**2
else
z(i)=sp01
end if
end do
case (12)
! rectangular emboss with specified slope
! sp01 is the z level
! sp02 and 03 are x1,x2
! sp04 and 05 are y1,y2
! sp06 is the thickness
! sp07 is the slope
m=tan(sp07*pi/180.d0)
x1a=sp02+(sp06/(2.d0*m))
x1b=sp02-(sp06/(2.d0*m))
x2a=sp03+(sp06/(2.d0*m))
x2b=sp03-(sp06/(2.d0*m))
y1a=sp04+(sp06/(2.d0*m))
y1b=sp04-(sp06/(2.d0*m))
y2a=sp05+(sp06/(2.d0*m))
y2b=sp05-(sp06/(2.d0*m))
do i=1,ns
if (x(i).le.x1a) then
zx=sp01
elseif (x(i).le.x1b) then
zx=sp01-(((x(i)-x1a)/(-sp06/m))*sp06)
elseif (x(i).le.x2a) then
zx=sp01-sp06
elseif (x(i).le.x2b) then
zx=sp01-((1.d0-((x(i)-x2a)/(-sp06/m)))*sp06)
!zx=sp01-((1.d0-((x(i)-x2a)/(sp06/m)))*sp06)
else
zx=sp01
endif
if (y(i).le.y1a) then
z(i)=sp01
elseif (y(i).le.y1b) then
z(i)=min(sp01-(((y(i)-y1a)/(-sp06/m))*sp06),zx)
elseif (y(i).le.y2a) then
z(i)=min(sp01-sp06,zx)
elseif (y(i).le.y2b) then
z(i)=min(sp01-((1.d0-((y(i)-y2a)/(-sp06/m)))*sp06),zx)
else
z(i)=sp01
endif
end do
case (13)
! a 3D embankment of finite length and with 2 kinks
! sp01 is the z level
! sp02 is the z level 2
! sp03 is slope angle (psi)
! sp04 is x value for the center of the slope (x0)
! sp05 is the x width of the narrow end of the plateau (y0w)
! sp06 is the y position for the start of the slope (y0)
! sp07 is the y position of the first kink (y1)
! sp08 is the angle of the first kink (theta1)
! sp09 is the y position of the second kink (y2)
! sp10 is the angle of the second kink (theta2)
! sp11 is the reference elevation (elevation of the lower outer dip change)
! sp12 is the dip angle of the lower outer dip panel
! sp13 is the reference elevation (elevation of the upper outer flat)
! sp14 is the dip angle of the upper outer dip panel
dz1=sp02-sp01
!dz2=sp01-sp11
dz2=sp01-sp11
dz3=sp01-sp13
dz4=dz3-dz2
psi=sp03*pi/180.d0
!psi2=sp12*pi/180.d0
psi2=sp12*pi/180.d0
psi3=sp14*pi/180.d0
theta1=sp08*pi/180.d0
theta2=sp10*pi/180.d0
psis=sign(tan(psi),dz1)
psi2s=sign(tan(psi2),dz2)
psi3s=sign(tan(psi3),dz3)
theta1s=tan(theta1)
theta2s=tan(theta2)
y0s=1.d0
y1s=tan(theta1/2.d0)
y2s=tan((theta1+theta2)/2.d0)
x0=sp04
y0w=sp05+dz1/(2.d0*psis)
y0=sp06+dz1/(2.d0*psis)
y0a=y0+dz1/(2.d0*psis)
y0b=y0-dz1/(2.d0*psis)
y1=sp07+(dz1/(2.d0*psis))*y1s
y2=sp09+(dz1/(2.d0*psis))*y2s
m1=((y2+(dz1/(2.d0*psis))*y1s)-(y1+(dz1/(2.d0*psis))*y1s))/&
((x0-y0w-(y1-y2)*theta1s-(dz1/(2.d0*psis)))-(x0-y0w-(dz1/(2.d0*psis))))
bb1=(y1+(dz1/(2.d0*psis))*y1s)-m1*(x0-y0w-(dz1/(2.d0*psis)))
a1=-m1
b1=1
c1=-bb1
m2=((y2-(dz1/(2.d0*psis))*y1s)-(y1-(dz1/(2.d0*psis))*y1s))/&
((x0-y0w-(y1-y2)*theta1s+(dz1/(2.d0*psis)))-(x0-y0w+(dz1/(2.d0*psis))))
bb2=(y1-(dz1/(2.d0*psis))*y1s)-m1*(x0-y0w+(dz1/(2.d0*psis)))
a2=-m2
b2=1
c2=-bb2
m3=((y2-(dz1/(2.d0*psis))*y1s)-(y1-(dz1/(2.d0*psis))*y1s))/&
((x0+y0w+(y1-y2)*theta1s-(dz1/(2.d0*psis)))-(x0+y0w-(dz1/(2.d0*psis))))
bb3=(y1-(dz1/(2.d0*psis))*y1s)-m3*(x0+y0w-(dz1/(2.d0*psis)))
a3=-m3
b3=1
c3=-bb3
m4=((y2+(dz1/(2.d0*psis))*y1s)-(y1+(dz1/(2.d0*psis))*y1s))/&
((x0+y0w+(y1-y2)*theta1s+(dz1/(2.d0*psis)))-(x0+y0w+(dz1/(2.d0*psis))))
bb4=(y1+(dz1/(2.d0*psis))*y1s)-m3*(x0+y0w+(dz1/(2.d0*psis)))
a4=-m4
b4=1
c4=-bb4
do i=1,ns
z(i)=sp01
if (y(i).ge.-y0s*x(i)+(y0+(x0-y0w)) .and. y(i).ge.y0s*x(i)+(y0-(x0+y0w))) then
if (y(i).ge.y0a) then
if (y(i).ge.y0a+(dz2/psi2s)) then
if (dz4.gt.eps) then
z(i)=max(sp01-dz2-psi3s*(y(i)-(y0a+(dz2/psi2s))),sp13)
else
z(i)=min(sp01-dz2-psi3s*(y(i)-(y0a+(dz2/psi2s))),sp13)
endif
else
z(i)=sp01-psi2s*(y(i)-y0a)
endif
! if (dz2.gt.eps) then
! z(i)=max(sp01-psi2s*(y(i)-y0a),sp11)
! else
! z(i)=min(sp01-psi2s*(y(i)-y0a),sp11)
! endif
elseif (y(i).ge.y0b) then
z(i)=sp01+psis*(y0a-y(i))
elseif (y(i).lt.y0b) then
z(i)=sp02
endif
elseif (y(i).ge.-y1s*x(i)+(y1s*(x0-y0w)+y1) .and. y(i).ge.y1s*x(i)+(-y1s*(x0+y0w)+y1)) then
if (x(i).lt.((x0-y0w)-(dz1/(2.d0*psis)))) then
if (x(i).lt.((x0-y0w)-(dz1/(2.d0*psis)))-(dz2/psi2s)) then
if (dz4.gt.eps) then
z(i)=max(sp01-dz2-psi3s*(((x0-y0w)-(dz1/(2.d0*psis))-(dz2/psi2s))-x(i)),sp13)
else
z(i)=min(sp01-dz2-psi3s*(((x0-y0w)-(dz1/(2.d0*psis))-(dz2/psi2s))-x(i)),sp13)
endif
else
z(i)=sp01-psi2s*(((x0-y0w)-(dz1/(2.d0*psis)))-x(i))
endif
! if ((dz2).gt.eps) then
! z(i)=max(sp01-psi2s*(((x0-y0w)-(dz1/(2.d0*psis)))-x(i)),sp11)
! else
! z(i)=min(sp01-psi2s*(((x0-y0w)-(dz1/(2.d0*psis)))-x(i)),sp11)
! endif
elseif (x(i).ge.((x0-y0w)-(dz1/(2.d0*psis))) .and. x(i).lt.((x0-y0w)+(dz1/(2.d0*psis)))) then
z(i)=sp01+psis*(x(i)-((x0-y0w)-(dz1/(2.d0*psis))))
elseif (x(i).ge.((x0-y0w)+(dz1/(2.d0*psis))) .and. x(i).lt.((x0+y0w)-(dz1/(2.d0*psis)))) then
z(i)=sp02
elseif (x(i).ge.((x0+y0w)-(dz1/(2.d0*psis))) .and. x(i).lt.((x0+y0w)+(dz1/(2.d0*psis)))) then
z(i)=sp01+psis*(x(i)-((x0+y0w)+(dz1/(2.d0*psis))))
else
if (x(i).lt.((x0+y0w)+(dz1/(2.d0*psis)))+(dz2/psi2s)) then
z(i)=sp01-psi2s*(x(i)-((x0+y0w)+(dz1/(2.d0*psis))))
else
if (dz4.gt.eps) then
z(i)=max(sp01-dz2-psi3s*(x(i)-((x0+y0w)+(dz1/(2.d0*psis))+(dz2/psi2s))),sp13)
else
z(i)=min(sp01-dz2-psi3s*(x(i)-((x0+y0w)+(dz1/(2.d0*psis))+(dz2/psi2s))),sp13)
endif
endif
! if ((dz2).gt.eps) then
! z(i)=max(sp01-psi2s*(x(i)-((x0+y0w)+(dz1/(2.d0*psis)))),sp11)
! else
! z(i)=min(sp01-psi2s*(x(i)-((x0+y0w)+(dz1/(2.d0*psis)))),sp11)
! endif
endif
elseif (y(i).ge.-y2s*x(i)+(y2s*(x0-y0w-(y1-y2)*theta1s)+y2) .and. y(i).ge.y2s*x(i)+(-y2s*(x0+y0w+(y1-y2)*theta1s)+y2)) then
if (y(i).ge.m1*x(i)+bb1) then
dist=abs(a1*x(i)+b1*y(i)+c1)/sqrt(a1**2.d0+b1**2.d0)
if (dist.lt.(dz2/psi2s)) then
z(i)=sp01-psi2s*dist
else
if (dz4.gt.eps) then
z(i)=max(sp01-dz2-psi3s*(dist-(dz2/psi2s)),sp13)
else
z(i)=min(sp01-dz2-psi3s*(dist-(dz2/psi2s)),sp13)
endif
endif
! if ((dz2).gt.eps) then
! z(i)=max(sp01-psi2s*dist,sp11)
! else
! z(i)=min(sp01-psi2s*dist,sp11)
! endif
elseif (y(i).ge.m2*x(i)+bb2 .and. y(i).lt.m1*x(i)+bb1) then
dist=abs(a1*x(i)+b1*y(i)+c1)/sqrt(a1**2.d0+b1**2.d0)
z(i)=sp01+dist*psis
elseif (y(i).lt.m2*x(i)+bb2 .and. y(i).lt.m3*x(i)+bb3) then
z(i)=sp02
elseif (y(i).ge.m3*x(i)+bb3 .and. y(i).lt.m4*x(i)+bb4) then
dist=abs(a4*x(i)+b4*y(i)+c4)/sqrt(a4**2.d0+b4**2.d0)
z(i)=sp01+dist*psis
else
dist=abs(a4*x(i)+b4*y(i)+c4)/sqrt(a4**2.d0+b4**2.d0)
if (dist.lt.(dz2/psi2s)) then
z(i)=sp01-psi2s*dist
else
if (dz4.gt.eps) then
z(i)=max(sp01-dz2-psi3s*(dist-(dz2/psi2s)),sp13)
else
z(i)=min(sp01-dz2-psi3s*(dist-(dz2/psi2s)),sp13)
endif
endif
! if ((dz2).gt.eps) then
! z(i)=max(sp01-psi2s*dist,sp11)
! else
! z(i)=min(sp01-psi2s*dist,sp11)
! endif
endif
elseif (x(i).lt.x0-y0w-(y1-y2)*theta1s-dz1/(2.d0*psis)) then
if (x(i).lt.x0-y0w-(y1-y2)*theta1s-dz1/(2.d0*psis)-(dz2/psi2s)) then
if (dz4.gt.eps) then
z(i)=max(sp01-dz2-psi3s*(((x0+y0w)+(y1-y2)*theta1s+(dz1/(2.d0*psis))+(dz2/psi2s))-x(i)),sp13)
else
z(i)=min(sp01-dz2-psi3s*(((x0+y0w)+(y1-y2)*theta1s+(dz1/(2.d0*psis))+(dz2/psi2s))-x(i)),sp13)
endif
else
z(i)=sp01-psi2s*((x0-y0w-(y1-y2)*theta1s-dz1/(2.d0*psis))-x(i))
endif
! if ((dz2).gt.eps) then
! z(i)=max(sp01-psi2s*((x0-y0w-(y1-y2)*theta1s-dz1/(2.d0*psis))-x(i)),sp11)
! else
! z(i)=min(sp01-psi2s*((x0-y0w-(y1-y2)*theta1s-dz1/(2.d0*psis))-x(i)),sp11)
! endif
elseif (x(i).ge.((x0-y0w-(y1-y2)*theta1s)-(dz1/(2.d0*psis))) .and. x(i).lt.((x0-y0w-(y1-y2)*theta1s)+(dz1/(2.d0*psis)))) then
z(i)=sp01+psis*(x(i)-((x0-y0w-(y1-y2)*theta1s)-(dz1/(2.d0*psis))))
elseif (x(i).ge.((x0-y0w-(y1-y2)*theta1s)+(abs(dz1)/(2.d0*psis))) .and. x(i).lt.((x0+y0w+(y1-y2)*theta1s)-(abs(dz1)/(2.d0*psis)))) then
z(i)=sp02
elseif (x(i).ge.((x0+y0w+(y1-y2)*theta1s)-(abs(dz1)/(2.d0*psis))) .and. x(i).lt.((x0+y0w+(y1-y2)*theta1s)+(abs(dz1)/(2.d0*psis)))) then
z(i)=sp01+psis*(x(i)-((x0+y0w+(y1-y2)*theta1s)+(dz1/(2.d0*psis))))
else
if (x(i).gt.x0+y0w+(y1-y2)*theta1s-dz1/(2.d0*psis)+(dz2/psi2s)) then
if (dz4.gt.eps) then
z(i)=max(sp01-dz2-psi3s*(x(i)-((x0+y0w)+(y1-y2)*theta1s+(dz1/(2.d0*psis))+(dz2/psi2s))),sp13)
else
z(i)=min(sp01-dz2-psi3s*(x(i)-((x0+y0w)+(y1-y2)*theta1s+(dz1/(2.d0*psis))+(dz2/psi2s))),sp13)
endif
else
z(i)=sp01-psi2s*(x(i)-(x0-y0w-(y1-y2)*theta1s-dz1/(2.d0*psis)))
endif
! if ((dz2).gt.eps) then
! z(i)=max(sp01-psi2s*(x(i)-(x0-y0w-(y1-y2)*theta1s-dz1/(2.d0*psis))),sp11)
! else
! z(i)=min(sp01-psi2s*(x(i)-(x0-y0w-(y1-y2)*theta1s-dz1/(2.d0*psis))),sp11)
! endif
endif
enddo
case default
call stop_run ('surface type not defined$')
end select
return
end