program perimetro_queima_1p0
  
  use portlib
  implicit none
    
    integer :: b    ! variable to open text files
    integer :: i    ! counting variable
    integer :: N    ! number of points of the star configuration
    real*8  :: w    ! web thickness [m]
    real*8  :: br   ! burning rate [m/s]
    real*8  :: r1   ! base radius [m]
    real*8  :: r2   ! tip radius [m]
    real*8  :: R    ! external radius [m]
    real*8  :: dt   ! time interval of computing S [s]
    real*8  :: etta ! angle [° and converted to rad]
    real*8  :: csi  ! secant fillet angle [° and converted to rad]
    real*8  :: pi   ! 3,141592...
    real*8  :: wx   ! web burned [m]
    real*8  :: S    ! perimeter of 1/N of the star [m]
    real*8  :: a    ! arc of the base [rad]
    real*8  :: Ystar ! auxiliary value to find the interface between zones 2 and 3 [m]
    
    ! the text file 'entries.out' is opened to the user write the proper inputs
    b = system('notepad inputs.out')
    
    ! reading the inputs
    open(1,file='inputs.out')
    read(1,*) N
    read(1,*) br
    read(1,*) w
    read(1,*) r1
    read(1,*) r2
    read(1,*) R
    read(1,*) etta
    read(1,*) csi
    read(1,*) dt
    close (1)
    
	pi = dacos(-1.0d0)
	csi = csi*pi/180.0d0
	etta = etta*pi/180.0d0

	! first some geometric requirements must be verified to guarantee the calculations
	if (csi > pi/N) then
        write(*,*) '---------'
		write(*,*) "Error: The angle csi can't be larger than pi/N!"
		write(*,*) '---------'
	elseif (etta < csi) then
        write(*,*) '---------'
		write(*,*) "Error: The angle etta can't be shorter than csi!"
	    write(*,*) '---------'
	else
		! two text files are opened to the user write the proper outputs
		! perimeter of burn function of the burned web
		open(2,file='outputs.out')
		! and perimeter of burn function of the time
		open(3,file='outputs2.out')
    
		a = pi/2.0d0 + csi - etta
		Ystar = (R-w-r1)*(dsin(csi)/dcos(etta)) - r1
		if (Ystar > w) then
		    write(*,*) '---------'
		    write(*,*) 'The web fraction is very small.'
		    write(*,*) 'Please increase external radius or web fraction'
		    write(*,*) '---------'
		else
		    wx = 0.0d0
		    S = (R-w+wx)*((pi/N)-csi) ! 1st part
		    S = S + (r1+wx)*a ! 2nd part
		    S = S + (R-w-r1)*(dsin(csi)/dsin(etta)) - (r1+r2)*dtan((pi/2.0d0) - etta) + (r2-wx)*((pi/2.0d0) - etta) ! 3rd part

		    write(2,*) wx, S*2.0d0*N
		    write(3,*) 0.0d0, S*2.0d0*N
        
		    do i = 1, 100000

			    wx = br*i*dt
            
			    ! zone 1 calculation
			    if (wx < r2) then
                
				    S = (R-w+wx)*((pi/N)-csi) ! 1st part
				    S = S + (r1+wx)*a ! 2nd part
				    S = S + (R-w-r1)*(dsin(csi)/dsin(etta)) - (r1+r2)*dtan((pi/2.0d0) - etta) + (r2-wx)*((pi/2.0d0) - etta) ! 3rd part
                
				    write(2,*) wx, S*2.0d0*N, '  z1'
				    write(3,*) i*dt, S*2.0d0*N
                
			    ! zone 2 calculation
			    elseif ((wx >= r2).and.(wx < Ystar)) then
                
				    S = (R-w+wx)*((pi/N)-csi) ! 1st part
				    S = S + (r1+wx)*a ! 2nd part
				    S = S + (R-w-r1)*(dsin(csi)/dsin(etta)) - (r1+wx)*dtan((pi/2.0d0) - etta) ! 3rd part
                
				    write(2,*) wx, S*2.0d0*N, '  z2'
				    write(3,*) i*dt, S*2.0d0*N
                
			    ! zone 3 calculation
			    elseif ((wx >= Ystar).and.(wx < w )) then
                
				    S = (R-w+wx)*((pi/N)-csi)
				    S = S + (r1+wx)*(csi + dasin( ( (R-w-r1)*dsin(csi) )  /(r1+wx))) ! 2nd part
                
				    write(2,*) wx, S*2.0d0*N, '  z3'
				    write(3,*) i*dt, S*2.0d0*N

			    ! zone 4 calculation
			    elseif (wx >= w) then
                
				    S = (r1+wx)*(csi + dasin(((R-w-r1)*dsin(csi))/(r1+wx)) - pi + dacos(((r1+wx)**2.0d0 + (R-r1-w)**2.0d0 - R*R)/(2*(r1+wx)*(R-r1-w))))
                
				    if(S < 0.d0) exit
                
				    write(2,*) wx, S*2.0d0*N, '  z4'
				    write(3,*) i*dt, S*2.0d0*N
                
			    end if    
            
		    end do
        
		    write(2,*) wx
        
		    close (2)
		    close (3)
        
		    !---------------------------------- 
		    open(9,file='sketch.txt') 
		    write(9,*) 'set grid'
		    write(9,*) "set xlabel 'burned web [m]'"
		    write(9,*) "set ylabel 'perimeter [m]'"
	        !write(9,*) "set size ratio -1"
	        !write(9,*) "set xrange[ ",0.0,": ",2.7d0*R," ]"
		    !write(9,*) "set yrange[ ",0.0,": ",15.0," ]"
        
		    write(9,*) "plot 'outputs.out' using 1:2 with lines"
		    !write(9,*) "plot 'outputs2.out' using 1:2 with lines"
		    close(9)
		    b = system('wgnuplot sketch.txt')
		    b = system('notepad outputs.out')

		    !---------------------------------- 
		    !open(9,file='sketch.txt') 
		    !write(9,*) 'set grid'
		    !write(9,*) "set xlabel 'time [s]'"
		    !write(9,*) "set ylabel 'burned web [m]'"
        
		    !write(9,*) "plot 'outputs2.out' using 1:2 with lines"
		    !close(9)
		    !b = system('wgnuplot sketch.txt')
		    !b = system('notepad outputs2.out')
        end if
	endif

end program perimetro_queima_1p0