Skip to content
Snippets Groups Projects
Commit a359efc7 authored by Dave Whipp's avatar Dave Whipp
Browse files

Fixed some typos in the new surface type 13

parent 5bae96cd
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 66 additions and 35 deletions
src/create_surfaces.f90
+ 66
− 35
View file @ a359efc7
......@@ -483,7 +483,8 @@ double precision sp01,sp02,sp03,sp04,sp05,sp06,sp07,sp08,sp09,sp10
integer i
double precision delta,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,c1,d1,e1,a2,b2,c2,d2,e2
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
!------------------------------------------------------------------------------|
......@@ -666,7 +667,8 @@ select case(surface_type)
! sp10 is the angle of the second kink (theta2)
psi=sp03*pi/180.d0
theta1=sp08*pi/180.d0
theta2=sp10*pi/180.d0
!theta2=sp10*pi/180.d0
theta2=0.d0
psis=tan(psi)
theta1s=tan(theta1)
theta2s=tan(theta2)
......@@ -676,53 +678,82 @@ select case(surface_type)
x0=sp04
y0w=sp05
y0=sp06
y0a=y0+(sp02/(2.d0*psis))
y0b=y0-(sp02/(2.d0*psis))
y0a=y0+(abs(sp02)/(2.d0*psis))
y0b=y0-(abs(sp02)/(2.d0*psis))
y1=sp07
y2=sp09
d1=((x0-y0w-(sp02/(2.d0*psis)))-(x0-y0w-(y1-y2)*theta1s-(sp02/(2.d0*psis))))
e1=((y1+(sp02/(2.d0*psis))*tan(theta1/2.d0))-(y2+(sp02/(2.d0*psis))*tan(theta1/2.d0)))
a1=-e1
b1=d1
c1=(-e1*(y1+(sp02/(2.d0*psis))*tan(theta1/2.d0))+d1*(x0-y0w-(sp02/(2.d0*psis))))
d2=((x0+y0w+(sp02/(2.d0*psis)))-(x0+y0w+(y1-y2)*theta1s+(sp02/(2.d0*psis))))
e2=((y1+(sp02/(2.d0*psis))*tan(theta1/2.d0))-(y2+(sp02/(2.d0*psis))*tan(theta1/2.d0)))
a2=-e2
b2=d2
c2=(e2*(y1+(sp02/(2.d0*psis))*tan(theta1/2.d0))-d2*(x0+y0w+(sp02/(2.d0*psis))))
m1=((y2+(abs(sp02)/(2.d0*psis))*tan(theta1/2.d0))-&
(y1+(abs(sp02)/(2.d0*psis))*tan(theta1/2.d0)))/&
((x0-y0w-(y1-y2)*theta1s-(abs(sp02)/(2.d0*psis)))-&
(x0-y0w-(abs(sp02)/(2.d0*psis))))
bb1=(y1+(abs(sp02)/(2.d0*psis))*tan(theta1/2.d0))-m1*(x0-y0w-(abs(sp02)/(2.d0*psis)))
a1=-m1
b1=1
c1=-bb1
m2=((y2-(abs(sp02)/(2.d0*psis))*tan(theta1/2.d0))-&
(y1-(abs(sp02)/(2.d0*psis))*tan(theta1/2.d0)))/&
((x0-y0w-(y1-y2)*theta1s+(abs(sp02)/(2.d0*psis)))-&
(x0-y0w+(abs(sp02)/(2.d0*psis))))
bb2=(y1-(abs(sp02)/(2.d0*psis))*tan(theta1/2.d0))-m1*(x0-y0w+(abs(sp02)/(2.d0*psis)))
a2=-m2
b2=1
c2=-bb2
m3=((y2-(abs(sp02)/(2.d0*psis))*tan(theta1/2.d0))-&
(y1-(abs(sp02)/(2.d0*psis))*tan(theta1/2.d0)))/&
((x0+y0w+(y1-y2)*theta1s-(abs(sp02)/(2.d0*psis)))-&
(x0+y0w-(abs(sp02)/(2.d0*psis))))
bb3=(y1-(abs(sp02)/(2.d0*psis))*tan(theta1/2.d0))-m3*(x0+y0w-(abs(sp02)/(2.d0*psis)))
a3=-m3
b3=1
c3=-bb3
m4=((y2+(abs(sp02)/(2.d0*psis))*tan(theta1/2.d0))-&
(y1+(abs(sp02)/(2.d0*psis))*tan(theta1/2.d0)))/&
((x0+y0w+(y1-y2)*theta1s+(abs(sp02)/(2.d0*psis)))-&
(x0+y0w+(abs(sp02)/(2.d0*psis))))
bb4=(y1+(abs(sp02)/(2.d0*psis))*tan(theta1/2.d0))-m3*(x0+y0w+(abs(sp02)/(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).lt.y0a .and. y(i).ge.y0b) then
z(i)=sp01+(((y(i)-y0a)/(-sp02/psis))*sp02)
z(i)=sp01+(((y(i)-y0a)/(-abs(sp02)/psis))*sp02)
elseif (y(i).lt.y0b) then
z(i)=sp01+sp02
z(i)=sp01+sp02+(x(i)-((x0-y0w-(y1-y2)*theta1s)+(abs(sp02)/(2.d0*psis))))*sp10
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).ge.((x0-y0w)-(sp02/(2.d0*psis))) .and. x(i).lt.((x0-y0w)+(sp02/(2.d0*psis)))) then
z(i)=sp01+(((x(i)-((x0-y0w)-(sp02/(2.d0*psis))))/(sp02/psis))*sp02)
elseif (x(i).ge.((x0-y0w)+(sp02/(2.d0*psis))) .and. x(i).lt.((x0+y0w)-(sp02/(2.d0*psis)))) then
z(i)=sp01+sp02
elseif (x(i).ge.((x0+y0w)-(sp02/(2.d0*psis))) .and. x(i).lt.((x0+y0w)+(sp02/(2.d0*psis)))) then
z(i)=sp01+(((x(i)-((x0+y0w)+(sp02/(2.d0*psis))))/(-sp02/psis))*sp02)
if (x(i).ge.((x0-y0w)-(abs(sp02)/(2.d0*psis))) .and. x(i).lt.((x0-y0w)+(abs(sp02)/(2.d0*psis)))) then
z(i)=sp01+(((x(i)-((x0-y0w)-(abs(sp02)/(2.d0*psis))))/(abs(sp02)/psis))*sp02)
elseif (x(i).ge.((x0-y0w)+(abs(sp02)/(2.d0*psis))) .and. x(i).lt.((x0+y0w)-(abs(sp02)/(2.d0*psis)))) then
z(i)=sp01+sp02+(x(i)-((x0-y0w-(y1-y2)*theta1s)+(abs(sp02)/(2.d0*psis))))*sp10
elseif (x(i).ge.((x0+y0w)-(abs(sp02)/(2.d0*psis))) .and. x(i).lt.((x0+y0w)+(abs(sp02)/(2.d0*psis)))) then
z(i)=sp01+(((x(i)-((x0+y0w)+(abs(sp02)/(2.d0*psis))))/(-abs(sp02)/psis))*sp02)
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.(theta1s*x(i)+(y1-theta1s*(x0-y0w)))-((1+y1s)*(sp02/(2.d0*psis))) .and. y(i).lt.(theta1s*x(i)+(y1-theta1s*(x0-y0w)))+((1+y1s)*(sp02/(2.d0*psis)))) then
!if (y(i).ge.(theta1s*x(i)+(y1-theta1s*(x0-y0w)))-((1+y1s)*(abs(sp02)/(2.d0*psis))) .and. y(i).lt.(theta1s*x(i)+(y1-theta1s*(x0-y0w)))+((1+y1s)*(abs(sp02)/(2.d0*psis)))) then
if (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/(sp02/psis)*sp02
elseif (y(i).lt.(theta1s*x(i)+(y1-theta1s*(x0-y0w)))-((1+y1s)*(sp02/(2.d0*psis))) .and. y(i).lt.(-theta1s*x(i)+(y1+theta1s*(x0+y0w)))-((1+y1s)*(sp02/(2.d0*psis)))) then
z(i)=sp01+sp02
elseif (y(i).ge.(-theta1s*x(i)+(y1+theta1s*(x0+y0w)))-((1+y1s)*(sp02/(2.d0*psis))) .and. y(i).lt.(-theta1s*x(i)+(y1+theta1s*(x0+y0w)))+((1+y1s)*(sp02/(2.d0*psis)))) then
dist=abs(a2*x(i)+b2*y(i)+c2)/sqrt(a2**2.d0+b2**2.d0)
z(i)=sp01+dist/(sp02/psis)*sp02
z(i)=sp01+dist/(abs(sp02)/psis)*sp02
!elseif (y(i).lt.(theta1s*x(i)+(y1-theta1s*(x0-y0w)))-((1+y1s)*(abs(sp02)/(2.d0*psis))) .and. y(i).lt.(-theta1s*x(i)+(y1+theta1s*(x0+y0w)))-((1+y1s)*(abs(sp02)/(2.d0*psis)))) then
elseif (y(i).lt.m2*x(i)+bb2 .and. y(i).lt.m3*x(i)+bb3) then
z(i)=sp01+sp02+(x(i)-((x0-y0w-(y1-y2)*theta1s)+(abs(sp02)/(2.d0*psis))))*sp10
!elseif (y(i).ge.(-theta1s*x(i)+(y1+theta1s*(x0+y0w)))-((1+y1s)*(abs(sp02)/(2.d0*psis))) .and. y(i).lt.(-theta1s*x(i)+(y1+theta1s*(x0+y0w)))+((1+y1s)*(abs(sp02)/(2.d0*psis)))) then
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/(abs(sp02)/psis)*sp02
endif
elseif (x(i).ge.((x0-y0w-(y1-y2)*theta1s)-(sp02/(2.d0*psis))) .and. x(i).lt.((x0-y0w-(y1-y2)*theta1s)+(sp02/(2.d0*psis)))) then
z(i)=sp01+(((x(i)-((x0-y0w-(y1-y2)*theta1s)-(sp02/(2.d0*psis))))/(sp02/psis))*sp02)
elseif (x(i).ge.((x0-y0w-(y1-y2)*theta1s)+(sp02/(2.d0*psis))) .and. x(i).lt.((x0+y0w+(y1-y2)*theta1s)-(sp02/(2.d0*psis)))) then
z(i)=sp01+sp02
elseif (x(i).ge.((x0+y0w+(y1-y2)*theta1s)-(sp02/(2.d0*psis))) .and. x(i).lt.((x0+y0w+(y1-y2)*theta1s)+(sp02/(2.d0*psis)))) then
z(i)=sp01+(((x(i)-((x0+y0w+(y1-y2)*theta1s)+(sp02/(2.d0*psis))))/(-sp02/psis))*sp02)
elseif (x(i).ge.((x0-y0w-(y1-y2)*theta1s)-(abs(sp02)/(2.d0*psis))) .and. x(i).lt.((x0-y0w-(y1-y2)*theta1s)+(abs(sp02)/(2.d0*psis)))) then
z(i)=sp01+(((x(i)-((x0-y0w-(y1-y2)*theta1s)-(abs(sp02)/(2.d0*psis))))/(abs(sp02)/psis))*sp02)
elseif (x(i).ge.((x0-y0w-(y1-y2)*theta1s)+(abs(sp02)/(2.d0*psis))) .and. x(i).lt.((x0+y0w+(y1-y2)*theta1s)-(abs(sp02)/(2.d0*psis)))) then
z(i)=sp01+sp02+(x(i)-((x0-y0w-(y1-y2)*theta1s)+(abs(sp02)/(2.d0*psis))))*sp10
elseif (x(i).ge.((x0+y0w+(y1-y2)*theta1s)-(abs(sp02)/(2.d0*psis))) .and. x(i).lt.((x0+y0w+(y1-y2)*theta1s)+(abs(sp02)/(2.d0*psis)))) then
z(i)=sp01+(((x(i)-((x0+y0w+(y1-y2)*theta1s)+(abs(sp02)/(2.d0*psis))))/(-abs(sp02)/psis))*sp02)
else
z(i)=sp01
endif
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment