program RAO_SIMPL
use portlib  
implicit none

real*8, dimension(:), allocatable :: Ee
real*8, dimension(:,:), allocatable :: ang_n, ang_s

real*8 :: temp

real*8 :: E, thetan, thetas, Rt

integer :: i, j
integer :: n, m, nx

i = system("Notepad dados.txt")
open(1,file='dados.txt')
    read(1,*) nx
    read(1,*) m
    read(1,*) E
    read(1,*) Rt
    Rt = Rt*0.5d0
close(1)

n = 35

allocate(Ee(n),ang_n(n,5),ang_s(n,5))

! Leitura das tabelas de angulo inicial
open(1,file='n_novo.txt')
read(1,*) 

do i=1, n
    read(1,*) temp, (ang_n(i,j),j=1,5)
enddo
close(1)

! Leitura das tabelas de angulo final
open(1,file='s_novo.txt')
read(1,*)
do i=1, n
    read(1,*) Ee(i), (ang_s(i,j),j=1,5)
enddo
close(1)


! Encontrar o valor de ang_n e ang_s para o caso escolhido:
call interp(E, Ee, n, m, ang_n, ang_s, thetan, thetas)

! Encontra a geometria sino simplificada de RAO:
call gera_RAO (m, Rt, E, thetan, thetas, nx)


contains
!------------------------------------------------------------
subroutine interp(E, Ee, n, m, ang_n, ang_s, thetan, thetas)
real*8, intent(in) :: E

integer, intent(in) :: n, m

real*8,intent(in), dimension(n) :: Ee
real*8, intent(in), dimension(n,5) :: ang_n, ang_s

real*8, intent(out) :: thetan, thetas

integer :: i, j, flag, pontos, inv
real*8 :: soma

real*8, dimension(n,3) :: dist
! Dist(:,1) = razao de expansao
! Dist(:,2) = ang_n
! Dist(:,3) = ang_s

! Declaracao
do i=1, n
    dist(i,1) = abs(E-Ee(i))
    dist(i,2) = ang_n(i,m)
    dist(i,3) = ang_s(i,m)
enddo

! Ordenacao por ordem de distancia para o valor de E
do j=1, 1000
    flag = 0
    do i=2, n
        if(dist(i,1) < dist(i-1,1))then
            temp = dist(i-1,1)
            dist(i-1,1) = dist(i,1)
            dist(i,1) = temp
            
            temp = dist(i-1,2)
            dist(i-1,2) = dist(i,2)
            dist(i,2) = temp
            
            temp = dist(i-1,3)
            dist(i-1,3) = dist(i,3)
            dist(i,3) = temp
            
            flag = 1
            
        endif
    enddo
    if(flag == 0) exit
enddo

if(dist(1,1) < 1E-10) then

    thetan = dist(1,2)
    thetas = dist(1,3)
else

! Ponderacao pela distancia, usando 2 pontos
pontos = 2
thetan = 0
thetas = 0
soma = 0
do i=1, pontos
    thetan = thetan + dist(i,2)/dist(i,1)
    thetas = thetas + dist(i,3)/dist(i,1)
    soma = soma + 1.d0/dist(i,1)
enddo
thetan = thetan/soma
thetas = thetas/soma

endif

end subroutine interp
!------------------------------------------------------------
subroutine gera_RAO(m, Rt, epsilon, thetan, thetas, nx)
use portlib
integer, intent(in) :: m, nx
real*8, intent(in) :: Rt, epsilon, thetan, thetas

real*8 :: R1,R2, Re, x0, tn, ts, ttn, tts
real*8 :: alpha15, leng, lengf, xn, yn, xe, ye, P, S, Q, T, Lf
real*8 :: pi=dacos(-1.d0)

real*8 :: dx, dx1, dx2, dx3, x, xp, y
integer :: passo1, passo2, passo3

select case (m)
    case (1)
        Lf = 0.6d0
    case (2)
        Lf = 0.7d0
    case (3)
        Lf = 0.8d0
    case (4)
        Lf = 0.9d0
    case (5)
        Lf = 1.0d0
end select

!  R1 = 1.5*Rt;
!  R2 = 0.382*Rt;
!  Re = sqrt(epsilon)*Rt;
!  x0 = -sqrt(2.0)*R1/2.0;
!  tn = thetan * M_PI / 180.0;
!  ts = thetas * M_PI / 180.0;
!  ttn = tan(tn);
!  tts = tan(ts);
!  alpha15 = 15.0 * M_PI / 180.0;
!  leng = Rt * (sqrt(epsilon)-1.0);
!  leng = leng + R1*((1.0/cos(alpha15))-1.0);
!  leng = leng / tan(alpha15);
!  lengf = (leng * Lf / 100.0);
!  xn = R2*cos(M_PI_2 - tn);
!  yn = Rt+R2-R2*sin(M_PI_2 - tn);
!  xe = lengf - xn;
!  ye = Re - yn;
!  P = (ye*ttn + ye*tts - 2*xe*tts*ttn) / (2*ye - xe*ttn - xe*tts);
!  S = (sqr(ye-P*xe)*(ttn-P)) / (xe*ttn-ye);
!  Q = -S / (2*(ttn-P));
!  T = Q*Q;

  R1 = 1.5d0*Rt;
  R2 = 0.382d0*Rt;
  Re = dsqrt(epsilon)*Rt;
  x0 = -dsqrt(2.0d0)*R1/2.0d0;
  tn = thetan * pi / 180.0d0;
  ts = thetas * pi / 180.0d0;
  ttn = dtan(tn);
  tts = dtan(ts);
  alpha15 = 15.0d0 * pi / 180.0d0;
  leng = Rt * (dsqrt(epsilon)-1.0d0);
  leng = leng + R1*((1.0d0/dcos(alpha15))-1.0d0);
  leng = leng / dtan(alpha15);
  lengf = leng * Lf ;
  xn = R2*dcos(pi/2 - tn);
  yn = Rt+R2-R2*dsin(pi/2 - tn);
  xe = lengf - xn;
  ye = Re - yn;
  P = (ye*ttn + ye*tts - 2*xe*tts*ttn) / (2*ye - xe*ttn - xe*tts);
  S = ((ye-P*xe)**2*(ttn-P)) / (xe*ttn-ye);
  Q = -S / (2.d0*(ttn-P));
  T = Q*Q;
  
    dx = (lengf-x0)/nx

    passo3 = (lengf-xn)/dx
    dx3 = (lengf-xn)/passo3

    passo2 = (xn-0)/dx
    dx2 = (xn- 0)/passo2

    passo1 = (0-x0)/dx
    dx1 = (0-x0)/passo1

    open(1,file='resultados.txt')
    
    write(1,*) "#Rt= ", Rt
    write(1,*) "#Epsilon= ", Epsilon
    write(1,*) "#thetan= ", thetan
    write(1,*) "#thetas= ", thetas    
    write(1,*) "#R1= ", R1
    write(1,*) "#R2= ", R2
    write(1,*) "#Re= ", Re
    write(1,*) "#x0= ", x0
    write(1,*) "#tn= ", tn
    write(1,*) "#ts= ", ts
    write(1,*) "#ttn= ", ttn
    write(1,*) "#tts= ", tts
    write(1,*) "#alpha15= ", alpha15
    write(1,*) "#leng= ", leng
    write(1,*) "#lengf= ", lengf
    write(1,*) "#xn= ", xn
    write(1,*) "#yn= ", yn
    write(1,*) "#xe= ", xe
    write(1,*) "#ye= ", ye
    write(1,*) "#P= ", P
    write(1,*) "#S= ", S
    write(1,*) "#Q= ", Q
    write(1,*) "#T= ", T
    write(1,*) "#"
    write(1,*) "#Lt= ", lengf-x0
    write(1,*) "#dx=", dx
    write(1,*) "#"
    write(1,*) "#passo1=", passo1
    write(1,*) "#dx1=", dx1
    write(1,*) "#"
    write(1,*) "#passo2=", passo2
    write(1,*) "#dx2=", dx2
    write(1,*) "#"
    write(1,*) "#passo3=", passo3
    write(1,*) "#dx3=", dx3

  ! primeiro arco de circulo
  !de x = x0 até x = 0
  x = x0
  !passo1 = (0-x0)/(lengf-xn)
  !dx = (0-x0)/passo1
  do while(x < 0.d0) 
  y = Rt + R1 - dsqrt(R1**2-x**2);
  
  write(1,*) x, y, -y
  x = x + dx1
  enddo

  ! segundo arco de circulo
  !de x = 0 até xn
  x = 0
 ! passo1 = (lengf-xn)/(xn-x)
!  dx = (xn-x)/passo1
  do while (x < xn)
    y = Rt + R2 - dsqrt(R2**2-x**2);
    write(1,*) x, y, -y
    x = x + dx2
  enddo
  
  
  
  ! parábola
  ! de x = xn até lengtf-xn
  xp = 0
  x = xn
  !dx3 = (lengf-xn)/passo3
  do while (x < lengf)
  ! xp = 0 até lengtf
!  ppasso = (lengf-xn) / ipasso;
    x = xp + xn;
    y = P*xp + Q + dsqrt(S*xp+T) + yn;

    write(1,*) x, y, -y
    xp = xp + dx3
    enddo
    
    
    close(1)
    
!    i= system("Notepad resultados.txt")
    
    open(1,file='plot.gnu')
        
write(1,*) "set term png nocrop enhanced size 1280,720"
write(1,*) "set out 'perfil_rao.png'"
write(1,*) "set style data lines"
write(1,*) "set grid"
write(1,*) "set xlabel 'x'"
write(1,*) "set ylabel 'y'"
write(1,*) "set time"
write(1,*) "set size ratio -1"
write(1,*) "set title 'Perfil RAO'"
write(1,*) "plot 'resultados.txt' using 1:2 t'',\"
write(1,*) " 'resultados.txt' using 1:3 t'',"
    
    close(1)
    i = system("Wgnuplot plot.gnu")
    i = system("start perfil_rao.png")

endsubroutine gera_RAO
!------------------------------------------------------------
end program RAO_SIMPL