Natrural Cubic Spline

!Este Programa fue realizado por Lillieskold, Thomas Enrique
!Instituto Universitario Aeronautico - 2022
!Docentes: Giovacchini, Juan Pablo; Weht, German
!Asignatura: Computacion e Introduccion al Calculo Numerico
!-------------------------------------------------------------------------------

!--------------------------DECLARACION DEL MODULO-------------------------------
MODULE Datos
!-----------------------------DATOS DE ENTRADA----------------------------------
    INTEGER,PARAMETER::N_Puntos=26
    INTEGER,PARAMETER::N_Nodos=N_Puntos-1

!------------------------VARIABLES DEL PROGRAMA---------------------------------

    REAL(8),DIMENSION(4*N_Nodos,4*N_Nodos)::Mat_Coef
    REAL(8),DIMENSION(4*N_Nodos)::Mat_Indep
    REAL(8),DIMENSION(4*N_Nodos)::Coeficiente
    REAL(8),DIMENSION(2,4)::Mat_Inf
    REAL(8),DIMENSION(2,4)::Mat_Sup
    REAL(8),DIMENSION(4,8)::Mat_Int
    REAL(8),DIMENSION(4)::Mat_Int_indep
    REAL(8),DIMENSION(N_Puntos)::Xpoint
    REAL(8),DIMENSION(N_Puntos)::Ypoint

!-------------------------------------------------------------------------------
!PROCEDO A DECLARAR VARIABLES PARA ELIMINACION GAUSSIANA

    REAL(8),DIMENSION(4*N_Nodos,4*N_Nodos)::A
    REAL(8),DIMENSION(4*N_Nodos)::B
    REAL(8),DIMENSION(4*N_Nodos)::X
    INTEGER N_Filas, N_Columnas

END MODULE

PROGRAM Spline_Cubico_Normal
    USE DATOS
    IMPLICIT NONE

    CALL Lectura_Datos_Par_Puntos
    CALL Burbuja_Mejorado(Xpoint,Ypoint)
    CALL Matrices_Externas
    CALL Matriz_Coef_Indep
    CALL Eliminacion_Gaussiana
    CALL Dots_Graph
    CALL Write_Outs

END PROGRAM

!SUBRUTINAS

!-------------------------------LECTURA DE DATOS--------------------------------
Subroutine Lectura_Datos_Par_Puntos
    USE Datos
    IMPLICIT NONE
    INTEGER I
        OPEN (2,File="Coordenadas_Intrados.txt")
        READ(2,*)
        READ(2,*)
        DO I=1,N_Puntos
            READ(2,*)Xpoint(i),Ypoint(i)
        END DO
End Subroutine
!-------------------------------------------------------------------------------

!--------------------------BURBUJA MEJORADO-------------------------------------
SUBROUTINE Burbuja_Mejorado(Vector_x1,Vector_y1)
    USE DATOS
    IMPLICIT NONE
REAL(8),DIMENSION(N_Puntos)::Vector_x1
REAL(8),DIMENSION(N_Puntos)::Vector_y1
REAL(8)Aux_x
REAL(8)Aux_y
LOGICAL ACCESO
INTEGER I,J

Aux_x=0.d0
Aux_y=0.d0

DO I=1,N_Puntos-1
    ACCESO=.False.
    DO J=1, N_Puntos-i
        IF(Vector_x1(j)>Vector_x1(j+1))Then
            !Pivoteo Vector x
            Aux_x=Vector_x1(j)
            Vector_x1(j)=Vector_x1(j+1)
            Vector_x1(j+1)=Aux_x

            !Pivoteo Vector y
            Aux_y=Vector_y1(j)
            Vector_y1(j)=Vector_y1(j+1)
            Vector_y1(j+1)=Aux_y

            Acceso=.True.
        END IF
    END DO
        If(.not.(Acceso)) EXIT
END DO
END SUBROUTINE Burbuja_Mejorado
!-------------------------------------------------------------------------------

!--------------------------MATRICES EXTERNAS------------------------------------
SUBROUTINE Matrices_Externas
    USE DATOS
    IMPLICIT NONE
    INTEGER J
    Mat_Inf=0.d0
    Mat_Sup=0.d0

    DO J=1,4
        Mat_Inf(1,J)=Xpoint(1)**(4-j)
        Mat_Sup(1,J)=Xpoint(N_Nodos+1)**(4-J)
    END DO
        Mat_Inf(2,1)=6.d0*Xpoint(1)
        Mat_Sup(2,1)=6.d0*Xpoint(N_Nodos+1)
        Mat_Inf(2,2)=2.d0
        Mat_Sup(2,2)=2.d0
END SUBROUTINE
!-------------------------------------------------------------------------------

!-------------------------MATRIZ_COEF_INDEP-------------------------------------
SUBROUTINE Matriz_Coef_Indep
    USE DATOS
    IMPLICIT NONE
    INTEGER I,M                !Variables de control para la matriz de coeficientes

    Mat_Coef=0.d0
    Mat_Indep=0.d0
    M=1

    Mat_Coef(1:2,1:4)=Mat_Inf
    Mat_Indep(1)=Ypoint(1)
    Mat_Coef(4*N_Nodos-1:4*N_Nodos,4*N_Nodos-3:4*N_Nodos)=Mat_Sup
    Mat_Indep(4*N_Nodos - 1)=Ypoint(N_Nodos+1)

    DO I=3, (4*N_Nodos)-2,4
        M=M+1
        CALL Matriz(Mat_Int,Mat_Int_indep,M)
        Mat_Indep(i:i+3)=Mat_int_indep
        Mat_Coef(i:i+3,i-2:i+5)=Mat_Int
    END DO
END SUBROUTINE
!-------------------------------------------------------------------------------

!------------------------------MATRIZ_INTERNA-----------------------------------
SUBROUTINE Matriz(Mat,Vec,M)
    USE DATOS
    IMPLICIT NONE
    REAL(8),DIMENSION(4,8)::Mat
    REAL(8),DIMENSION(4)::Vec
    INTEGER J,M

    Mat=0.d0
    Vec=0.d0

        DO J=1,4
            Mat(1,J)=Xpoint(M)**(4-J)
        END DO

        DO J=5,8
            Mat(2,J)=Xpoint(M)**(8-J)
        END DO

            Mat(3,1)=3*Xpoint(M)**2
            Mat(3,2)=2*Xpoint(M)
            Mat(3,3)=1.d0
            Mat(3,5)= -Mat_Int(3,1)
            Mat(3,6)= -Mat_Int(3,2)
            Mat(3,7)= -Mat_Int(3,3)
            Mat(4,1)=6*Xpoint(M)
            Mat(4,2)=2.d0
            Mat(4,5)= -Mat_Int(4,1)
            Mat(4,6)= -Mat_Int(4,2)

            Vec(1)=Ypoint(m)
            Vec(2)=Ypoint(m)

END SUBROUTINE
!-------------------------------------------------------------------------------

!------------------------ELIMINACION_GAUSSIANA----------------------------------
SUBROUTINE Eliminacion_Gaussiana
    USE DATOS
    IMPLICIT NONE

    !PROCEDO A ASIGNAR VALORES A LAS VARIABLES PARA ELIMINACION GAUSSIANA

    A=Mat_Coef
    B=Mat_Indep
    X=Coeficiente
    N_Filas=4*N_Nodos
    N_Columnas=4*N_Nodos


    CALL Triangular_sup
    CALL Retrosustitucion


    Coeficiente=X

END SUBROUTINE

SUBROUTINE Triangular_sup
    USE DATOS
    IMPLICIT NONE
    REAL(8)Multiplicador
    INTEGER i,j,k
    DO i= 1 , N_Filas-1
        CALL PIVOTEO (i)
        DO j= i+1 , N_Filas
            Multiplicador= A(j,i)/A(i,i)
            A(j,i)=0                        !Con esto me aseguro que las componentes a la izquierda de la diagonal sean cero
            B(j)=B(j)-Multiplicador*B(i)
            DO K= i+1 , N_Filas
                A(j,k)= A(j,k)-Multiplicador*A(i,k)
            END DO
        END DO
    END DO
END SUBROUTINE Triangular_sup

SUBROUTINE Pivoteo(Ref_Columna)
    USE DATOS
    IMPLICIT NONE
REAL(8),DIMENSION(N_Columnas)::Aux
REAL(8) Ref_Coef, Aux_B
INTEGER Ref_Fila,Ref_Columna
INTEGER j,k

    Ref_Coef=A(Ref_Columna,Ref_Columna)                     !Coeficiente de referencia -> DIAGONAL PRINCIPAL
    Ref_Fila=Ref_Columna                                    !Variable que guarda la fila del coeficiente de referencia
    DO j= Ref_Columna+1,N_Filas                               !Recorre el subvector -> El valor es i+1 debido a que me interesan los valores por debajo de la diagonal principal
        IF (ABS(Ref_Coef)<ABS(A(j,Ref_Columna)))THEN
            Ref_Coef=A(j,Ref_Columna)                      !Hago que el valor de referencia para la comparacion sea el mayor
            Ref_Fila=j                                     !Fila en la cual esta el elemento de referencia
        END IF
    END DO                      !Cuando termina este bucle tengo en Ref_Fila la fila en la que se encuentra el valor maximo

    DO k=1,N_Columnas
    Aux(k)= A(Ref_Columna,k)            !Guardo la fila de referencia inicial completa - !Uso Ref_Columna debido a que es igual a la fila que me interesa
    END DO

    DO k=1,N_Columnas
    A(Ref_Columna,k)=A(Ref_Fila,k)           !Reemplazo la fila de referencia inicial por la que tiene el coeficiente de valor mas elevado
    END DO

    DO k=1,N_Columnas
    A(Ref_Fila,k)=Aux(k)      !Reescribo la fila de referencia inicial en la posicion de la que tenia el coeficiente de valor mas elevado
    END DO

    Aux_B=B(Ref_Fila)
    B(Ref_Fila)=B(Ref_Columna)
    B(Ref_Columna)=Aux_B
END SUBROUTINE Pivoteo

SUBROUTINE Retrosustitucion
    USE DATOS
    IMPLICIT NONE
    REAL(8)Producto
    INTEGER I,J
    DO i= N_Filas , 1 , -1
    Producto=0.d0
    DO j= i+1 , N_Filas
        Producto= Producto + (A(i,j)*X(j))
    END DO
    X(i)=(B(i)-Producto)/A(i,i)
END DO
END SUBROUTINE
!-------------------------------------------------------------------------------

!-------------------------------DOTS GRAPH--------------------------------------
Subroutine Dots_Graph
    USE Datos
    IMPLICIT NONE
    REAL(8) Func
    REAL(8)X_Point
    REAL(8)dx_Point
    INTEGER N_PointToPlot, Point_I

    N_PointToPlot=100
    dx_Point=(Maxval(Xpoint)-Minval(Xpoint))/(N_PointToPlot-1)
    Do Point_I=1,N_PointToPlot
        X_point= Minval(XPoint)+(Point_I-1)*dx_Point
    WRITE(50,*)X_Point,Func(X_point)
    End do
End Subroutine
!-------------------------------------------------------------------------------

REAL(8)FUNCTION Func(X_int)
    USE DATOS
    IMPLICIT NONE
    REAL(8)X_int
    REAL(8)a_i,b_i,c_i,d_i
    INTEGER I

    DO I=1,N_Nodos
        IF(X_int>Xpoint(i) .and. X_Int<Xpoint(i+1))THEN
            !Estoy en P_i
            a_i=Coeficiente(4*i-3)
            b_i=Coeficiente(4*i-2)
            c_i=Coeficiente(4*i-1)
            d_i=Coeficiente(4*i)
        END IF
        FUNC=a_i*X_int**3+b_i*X_int**2+c_i*X_int+d_i
    END DO
END FUNCTION


!--------------------------------WRITE OUTS-------------------------------------
SUBROUTINE Write_Outs
    USE DATOS
    IMPLICIT NONE
    INTEGER I,K,J
    OPEN(3,File="Write_Outs_Ejercicio_1_Final.txt")
    k=1
    J=0
        DO I=1,4*N_Nodos,4
            J=J+1
                Print*,"Coeficientes de P_",J, "[",(Coeficiente(I+K-1),K=1,4),"]"
                PRINT*, "--------------------------------------------------------------"
                WRITE(3,*) "Coeficientes de P_",J, "[",(Coeficiente(I+K-1),K=1,4),"]"
                WRITE(3,*)"-------------------------------------------------------------"
        END DO
END SUBROUTINE