commit 5f8f33347ffff4b1266a295b7757251baaf7cb88
parent eb6195eb78ec09de2749fbe08a2c69644a9ac0ca
Author: Julio C. Ramirez-Ceballos <julio@jcrcx.xyz>
Date: Wed, 4 Dec 2024 22:16:56 +0100
Added function to write output file. (and made some updates on the other stuff)
Diffstat:
| M | README | | | 16 | ++++++++++------ |
| M | clt.c | | | 49 | ++++++++++++++++++++++++++++--------------------- |
| M | clt.h | | | 3 | ++- |
| M | main.c | | | 77 | ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++----------- |
| M | rw.c | | | 14 | +++++++++----- |
| A | utils.c | | | 190 | +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
| M | utils.h | | | 279 | ++++++------------------------------------------------------------------------- |
7 files changed, 323 insertions(+), 305 deletions(-)
diff --git a/README b/README
@@ -13,14 +13,18 @@ optimization of the laminate.
The program uses Classical Laminate Theory (CLT) to compute the strain and
stress on each node of a given 2D mesh with some albitrary boundary conditions.
-**Input**:
+INPUT
+-----
+
+<!-- ARDO will use Leap motion for hand tracking -->
- Mesh generated from Gmsh. It must include defined boundary conditions.
- Materials. These are given on the `materials.dat` file
- Laminate. These are given on the `laminate.dat` file.
- - [ ] The laminate can be optimized later with **ARDO** using *kinematic
- draping theory*
-- Boundary conditions. These are given on the `bc.dat` file. This file links the
-defined physical objects on the mesh file to apply displacements and loads
-(mechanical or hygro-thermal).
+ - [ ] The laminate can be optimized later with **ARDO** using *kinematic
+ draping theory*
+ - [ ] Generate laminate as chromosomes
+- Boundary conditions. These are given on the `loads.dat` file. This file links
+ the defined physical objects on the mesh file to apply displacements and loads
+ (mechanical or hygro-thermal).
diff --git a/clt.c b/clt.c
@@ -6,7 +6,7 @@
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
-#include "utils.h"
+#include "utils.c"
#include "rw.c"
#include "clt.h"
@@ -419,9 +419,11 @@ double E_x(double ABD[6][6], double h) {
return Ex;
}
-double E_y(double ABD[6][6], double h) {
+void
+E_y(double ABD[6][6], double h, double Ey)
+{
double matrix1[6][6], matrix2[5][5];
- double det1, det2, Ey;
+ double det1, det2;
// Copy relevant parts of ABD to matrix1 for the first determinant
for (i = 0; i < 6; i++) {
@@ -433,7 +435,7 @@ double E_y(double ABD[6][6], double h) {
// Perform LU decomposition for matrix1
if (LUDecompose(&matrix1[0][0], 6, 6) < 0) {
fprintf(stderr, "Error during LU decomposition for matrix1.\n");
- return -1;
+ return;
}
// Calculate determinant for matrix1
@@ -468,7 +470,7 @@ double E_y(double ABD[6][6], double h) {
// Perform LU decomposition for matrix2
if (LUDecompose(&matrix2[0][0], 5, 5) < 0) {
fprintf(stderr, "Error during LU decomposition for matrix2.\n");
- return -1;
+ return;
}
// Calculate determinant for matrix2
@@ -477,15 +479,16 @@ double E_y(double ABD[6][6], double h) {
// Check for division by zero
if (det2 == 0) {
fprintf(stderr, "Error: Division by zero in E_x calculation, determinant of the submatrix is zero.\n");
- return -1;
+ return;
}
// Compute E_x
Ey = (det1 / det2) / h;
- return Ey;
}
-double G_xy(double ABD[6][6], double h) {
+double
+G_xy(double ABD[6][6], double h)
+{
double matrix1[6][6], matrix2[5][5];
double det1, det2, Gxy;
@@ -555,7 +558,9 @@ double G_xy(double ABD[6][6], double h) {
return Gxy;
}
-double nu_xy(double ABD[6][6], double h) {
+double
+nu_xy(double ABD[6][6], double h)
+{
double matrix1[5][5], matrix2[5][5];
double det1, det2, nxy;
@@ -640,7 +645,9 @@ double nu_xy(double ABD[6][6], double h) {
return nxy;
}
-double nu_yx(double ABD[6][6], double h) {
+double
+nu_yx(double ABD[6][6], double h)
+{
double matrix1[5][5], matrix2[5][5];
double det1, det2, nyx;
@@ -716,7 +723,7 @@ double nu_yx(double ABD[6][6], double h) {
// Check for division by zero
if (det2 == 0) {
- fprintf(stderr, "Error: Division by zero in E_y calculation, determinant of the submatrix is zero.\n");
+ fprintf(stderr, "Error: Division by zero in n_yx calculation, determinant of the submatrix is zero.\n");
return -1;
}
@@ -820,9 +827,9 @@ print_lcs_mcs(int num_layers,
double MS_stress_inf[3][MAX_LAYERS])
{
printf("\n");
- printf("+-------------------------+\n");
- printf("| Laminate System Strains |\n");
- printf("+-------------------------+\n");
+ printf("+------------------------------------+\n");
+ printf("| Laminate Coordinate System Strains |\n");
+ printf("+------------------------------------+\n");
printf("\n");
printf("Layer | Inferior (Z-) | Superior (Z+)\n");
printf("------|--------------------------------|----------------------------\n");
@@ -840,9 +847,9 @@ print_lcs_mcs(int num_layers,
}
printf("\n");
- printf("+-------------------------+\n");
- printf("| Material System Strains |\n");
- printf("+-------------------------+\n");
+ printf("+------------------------------------+\n");
+ printf("| Material Coordinate System Strains |\n");
+ printf("+------------------------------------+\n");
printf("\n");
printf("Layer | Inferior (Z-) | Superior (Z+)\n");
printf("------|--------------------------------|----------------------------\n");
@@ -860,9 +867,9 @@ print_lcs_mcs(int num_layers,
}
printf("\n");
- printf("+------------------------+\n");
- printf("| Material System Stress |\n");
- printf("+------------------------+\n");
+ printf("+-----------------------------------+\n");
+ printf("| Material Coordinate System Stress |\n");
+ printf("+-----------------------------------+\n");
printf("\n");
printf("Layer | Inferior (Z-) | Superior (Z+)\n");
printf("------|--------------------------------|----------------------------\n");
@@ -884,7 +891,7 @@ void
print_equivalent_properties(double ABD[6][6], double h)
{
double Ex = E_x(ABD, h);
- double Ey = E_y(ABD, h);
+ double Ey = 0; E_y(ABD, h, Ey);
double Gxy = G_xy(ABD, h);
double nxy = nu_xy(ABD, h);
double nyx = nu_yx(ABD, h);
diff --git a/clt.h b/clt.h
@@ -52,7 +52,7 @@ void maxstress_2D(MaterialProperties *mat_list, Laminate *lam,
double **stress_inf, double **stress_sup,
FSResult *fs_inf, FSResult *fs_sup);
double E_x(double ABD[6][6], double h);
-double E_y(double ABD[6][6], double h);
+void E_y(double ABD[6][6], double h, double Ey);
double G_xy(double ABD[6][6], double h);
double nu_xy(double ABD[6][6], double h);
double nu_yx(double ABD[6][6], double h);
@@ -64,3 +64,4 @@ void print_lcs_mcs(int num_layers,
double MS_stress_sup[3][MAX_LAYERS],
double MS_stress_inf[3][MAX_LAYERS]);
void print_equivalent_properties(double ABD[6][6], double h);
+/* ADD FUNCTION TO COMPUTE THE PRINCIPAL STRESS/STRAIN FOR EACH PLY */
diff --git a/main.c b/main.c
@@ -1,7 +1,4 @@
-/*
- * main.c
- *
- * Author: Julio C. Ramirez-Ceballos
+/* Author: Julio C. Ramirez-Ceballos
* Licence: ISC
* Description:
* This code is used to calculate the equivalent properties of a laminate
@@ -23,24 +20,82 @@
*/
#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
#include "clt.c"
-int main()
+void write_output
+(const char *filename, const char *content)
+{
+ FILE *file = fopen(filename, "w");
+ if (file == NULL) {
+ printf("Error opening file!\n");
+ exit(1);
+ }
+ fprintf(file, "%s", content);
+ fclose(file);
+}
+
+void print_usage() {
+ printf("Usage: ./main [-p | -w]\n");
+ printf(" -p : Print output to console\n");
+ printf(" -w : Write output to output.res file\n");
+}
+
+
+int main
+(int argc, char *argv[])
{
+ if (argc != 2 || (strcmp(argv[1], "-p") != 0 && strcmp(argv[1], "-w") != 0)) {
+ print_usage();
+ return 1;
+ }
+
+ int save_to_file = (strcmp(argv[1], "-w") == 0);
+
+ /* Improve function for saving results and debugging */
+ /* Option to save or print output */
+ char output[10000]; /* Assuming the output won't exceed this size */
+ int output_length = 0;
+
+ /* Redirect stdout to our buffer if we're saving to file */
+ if (save_to_file) {
+ freopen("output.res", "w", stdout);
+ }
+
int num_materials;
MaterialProperties *materials =
- read_material_properties("materials.dat", &num_materials, true);
+ read_material_properties("materials.dat", &num_materials, 1);
Laminate lam;
- read_laminate_data("laminate.dat", &lam, true);
+ read_laminate_data("laminate.dat", &lam, 1);
- double F[N] = {0}; // Initialize F with zeros
+ double F[6];
read_force_vector(F, "loads.dat");
- /* print, ABD, Z */
- clt(&lam, materials, num_materials, F, 1, 0 , 0);
+ /* Perform CLT calculation */
+ clt(&lam, materials, num_materials, F, 1, 0, 0);
+
+ /* If we're saving to file, we need to capture the output */
+ if (save_to_file) {
+ fclose(stdout); /* Close the file stream */
+ FILE *file = fopen("output.res", "r");
+ if (file == NULL) {
+ printf("Error reading file!\n");
+ exit(1);
+ }
+ while
+ (fgets(output + output_length,
+ sizeof(output) - output_length,
+ file) != NULL)
+ {
+ output_length += strlen(output + output_length);
+ }
+ fclose(file);
+ write_output("output.res", output);
+ }
- // Free allocated memory
+ /* Free allocated memory */
free(lam.thk);
free(lam.ang);
free(lam.mat_id);
diff --git a/rw.c b/rw.c
@@ -26,10 +26,13 @@ typedef struct {
int num_layers;
} Laminate;
-void read_laminate_data(const char *filename, Laminate *lam, bool print_details) {
+void read_laminate_data
+(const char *filename, Laminate *lam, bool print_details)
+{
FILE *file;
char line[100]; // Assuming each line won't exceed 100 characters
int layer_count = 0;
+ double h = 0;
// Attempt to open the file
file = fopen(filename, "r");
@@ -81,6 +84,8 @@ void read_laminate_data(const char *filename, Laminate *lam, bool print_details)
}
}
}
+ for(int i = 0; i < lam->num_layers; i++) h += lam->thk[i];
+
// Print the laminate details if requested
if (print_details) {
@@ -89,6 +94,7 @@ void read_laminate_data(const char *filename, Laminate *lam, bool print_details)
printf("| Laminate Configuration |\n");
printf("+------------------------+\n\n");
printf("Number of Layers: %d\n", lam->num_layers);
+ printf("Laminate total thickness = %.4e\n", h);
printf("Layer | Thickness (mm) | Angle (°) | Material ID\n");
printf("------|----------------|-----------|------------\n");
for (int i = 0; i < lam->num_layers; i++) {
@@ -103,10 +109,8 @@ void read_laminate_data(const char *filename, Laminate *lam, bool print_details)
fclose(file);
}
-MaterialProperties*
-read_material_properties(const char *filename,
- int *n_mat,
- bool print_details)
+MaterialProperties *read_material_properties
+(const char *filename, int *n_mat, bool print_details)
{
FILE *file;
char line[MAX_LINE_LENGTH];
diff --git a/utils.c b/utils.c
@@ -0,0 +1,190 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#include "utils.h"
+
+/* Factorial for Taylor series */
+static double factorial(int n) {
+ if (n == 0 || n == 1) return 1;
+ return n * factorial(n - 1);
+}
+
+/* Sine function using Taylor series */
+double Sin
+(double x)
+{
+ double sum = 0;
+ int n;
+ x = fmod(x, 2 * PI); /* Reduce x to [0, 2π) */
+ for (n = 0; n < 10; ++n) {
+ sum += pow(-1, n) * pow(x, 2*n+1) / factorial(2*n+1);
+ }
+ return sum;
+}
+
+/* Cosine function using Taylor series */
+double
+Cos(double x) {
+ double sum = 0;
+ int n;
+ x = fmod(x, 2 * PI); /* Reduce x to [0, 2π) */
+ for (n = 0; n < 10; ++n) {
+ sum += pow(-1, n) * pow(x, 2*n) / factorial(2*n);
+ }
+ return sum;
+}
+
+/* Matrix multiplication */
+void
+matmul(double A[3][3], double B[3][3], double result[3][3])
+{
+ for (int i = 0; i < 3; i++) {
+ for (int j = 0; j < 3; j++) {
+ result[i][j] = 0;
+ for (int k = 0; k < 3; k++) {
+ result[i][j] += A[i][k] * B[k][j];
+ }
+ }
+ }
+}
+
+/* Vector-matrix multiplication */
+void
+dot_product(double result[3], double A[3][3], double B[3])
+{
+ for(int i = 0; i < 3; i++) {
+ result[i] = 0;
+ for(int j = 0; j < 3; j++) {
+ result[i] += A[i][j] * B[j];
+ }
+ }
+}
+
+/* Matrix transposition */
+void
+transpose(double A[3][3], double Ai[3][3])
+{
+ for (int j = 0; j < 3; j++) {
+ for (int k = 0; k < 3; k++) {
+ Ai[j][k] = A[k][j];
+ }
+ }
+}
+
+/* Display matrix in row or column major order */
+void
+dispmat(double *A, int rows, int cols, bool transpose)
+{
+ if (transpose) {
+ for (int j = 0; j < cols; j++) {
+ printf("%4d ", j + 1);
+ for (int i = 0; i < rows; i++) {
+ printf("%10.2e", *(A + i * cols + j));
+ }
+ printf("\n");
+ }
+ } else {
+ for (int i = 0; i < rows; i++) {
+ for (int j = 0; j < cols; j++) {
+ printf("%10.2e", *(A + i * cols + j));
+ }
+ printf("\n");
+ }
+ }
+}
+
+/* Static variables for pivot */
+static int *pivot = NULL;
+
+/* LU Decomposition */
+int LUDecompose(double *amat, int n, int numcols) {
+ if (pivot != NULL) free(pivot);
+ pivot = malloc(n * sizeof(int));
+ if (!pivot) {
+ fprintf(stderr, "Error in LUDecompose - malloc\n");
+ return -2;
+ }
+
+ for (int i = 0; i < n; i++) {
+ pivot[i] = i;
+ }
+
+ int sign = 1;
+ double *dptr1, *dptr2, dtmp1;
+
+ for (int i = 0; i < n-1; i++) {
+ dptr1 = amat + i * numcols + i;
+ int k = i;
+ for (int j = i + 1; j < n; j++) {
+ dptr2 = amat + j * numcols + i;
+ if (fabs(*dptr2) > fabs(*dptr1)) {
+ dptr1 = dptr2;
+ k = j;
+ }
+ }
+
+ if (fabs(*dptr1) < SMALLEST_PIVOT) {
+ fprintf(stderr, "Error in LUDecompose - matrix is singular\n");
+ return -1;
+ }
+
+ if (k != i) {
+ int temp = pivot[i]; pivot[i] = pivot[k]; pivot[k] = temp;
+ sign *= -1;
+ for (int j = 0; j < n; j++) {
+ dtmp1 = amat[i * numcols + j];
+ amat[i * numcols + j] = amat[k * numcols + j];
+ amat[k * numcols + j] = dtmp1;
+ }
+ }
+
+ for (int j = i+1; j < n; j++) {
+ dtmp1 = amat[j * numcols + i] / *dptr1;
+ amat[j * numcols + i] = dtmp1;
+ for (int k = i+1; k < n; k++) {
+ amat[j * numcols + k] -= dtmp1 * amat[i * numcols + k];
+ }
+ }
+ }
+ /*pivot[n-1] = sign;*/
+ *(pivot+n-1) = sign;
+ return 1;
+}
+
+/* Solve system with LU decomposition */
+void LUSolve(double *amat, double *b, int n, int numcols) {
+ double dtmp1;
+ int i, j;
+
+ for (i = 0; i < n; i++) {
+ j = pivot[i];
+ if (j != i) {
+ dtmp1 = b[i]; b[i] = b[j]; b[j] = dtmp1;
+ }
+ }
+
+ for (i = 0; i < n; i++) {
+ dtmp1 = b[i];
+ for (j = 0; j < i; j++) {
+ dtmp1 -= amat[i * numcols + j] * b[j];
+ }
+ b[i] = dtmp1;
+ }
+
+ for (i = n-1; i >= 0; i--) {
+ dtmp1 = b[i];
+ for (j = n-1; j > i; j--) {
+ dtmp1 -= amat[i * numcols + j] * b[j];
+ }
+ b[i] = dtmp1 / amat[i * numcols + i];
+ }
+}
+
+/* Compute determinant from LU decomposition */
+double LUDeterminant(double *amat, int n, int numcols) {
+ double det = 1.0;
+ for (int i = 0; i < n; i++) {
+ det *= amat[i * (numcols + 1)];
+ }
+ return det * pivot[n - 1];
+}
diff --git a/utils.h b/utils.h
@@ -1,269 +1,26 @@
-/*
- * utils.h
- */
+#ifndef UTILS_H
+#define UTILS_H
#include <math.h>
#include <stdbool.h>
#define PI 3.14159265358979323846
-#define DEGTORAD 0.017453292519943295 /* pi/180 */
-#define CONST13 0.333333333333333
+#define DEGTORAD (PI/180)
+#define CONST13 (1.0/3.0)
#define MAT_ELEM(mat, row, col) mat[row][col]
-#define SMALLEST_PIVOT 1.0e-5 /* Smallest allowed non-singular pivot. */
+#define SMALLEST_PIVOT (1e-5)
#define SIZE 6
-static int *pivot = NULL; /* Pivot vector. */
-
-/* Simple factorial function for Taylor series */
-double factorial(int n) {
- if (n == 0 || n == 1) return 1;
- return n * factorial(n - 1);
-}
-
-/* Sine function using Taylor series */
-double Sin(double x) {
- double term, sum = 0;
- int n;
- x = fmod(x, 2 * PI); /* Reduce x to be within[0, 2π) */
-
- for (n = 0; n < 10; ++n) { /* Adjust the number of terms for precision */
- term = pow(-1, n) * pow(x, 2*n+1) / factorial(2*n+1);
- sum += term;
- }
- return sum;
-}
-
-/* Cosine function using Taylor series */
-double Cos(double x) {
- double term, sum = 0;
- int n;
- x = fmod(x, 2 * PI); /* Reduce x to be within[0, 2π) */
-
- /* Adjust the number of terms for precision */
- for (n = 0; n < 10; ++n) {
- term = pow(-1, n) * pow(x, 2*n) / factorial(2*n);
- sum += term;
- }
- return sum;
-}
-
-void matmul(double A[3][3], double B[3][3], double result[3][3])
-{
- for (int i = 0; i < 3; i++) {
- for (int j = 0; j < 3; j++) {
- result[i][j] = 0.0;
- for (int k = 0; k < 3; k++) {
- result[i][j] += A[i][k] * B[k][j];
- }
- }
- }
-}
-
-/*CORRIGER*/
-void dot_product(double result[3], double A[3][3], double B[3])
-{
- for(int i = 0; i < 3; i++) {
- result[i] = 0;
- for(int j = 0; j < 3; j++) {
- result[i] += A[i][j] * B[j];
- }
- }
-}
-
-void transpose(double A[3][3], double Ai[3][3])
-{
- for (int j = 0; j < 3; j++) {
- for (int k = 0; k < 3; k++) {
- Ai[j][k] = A[k][j];
- }
- }
-}
-
-void dispmat(double *A, int rows, int cols, bool transpose)
-{
- if (transpose) {
- for (int j = 0; j < cols; j++) {
- printf("%4d ", j + 1);
- for (int i = 0; i < rows; i++) {
- printf("%10.2e", *(A + i * cols + j));
- }
- printf("\n");
- }
- } else {
- for (int i = 0; i < rows; i++) {
- for (int j = 0; j < cols; j++) {
- printf("%10.2e", *(A + i * cols + j));
- }
- printf("\n");
- }
- }
-
-}
-
-/*
- * This file contains routines for solving an nxn system of linear equations
- * using LU decomposition. The implementation is based on material in
- * "Numerical Methods - A Software Approach" by Johnston.
- * Function LUDecompose performs the decomposition step, LUSolve solves a
- * system given a particular right hand side vector b, and LUDetermiant will
- * return the determinant of a matrix which has been decomposed to LU form.
- */
-
-int LUDecompose(double *amat, int n, int numcols)
-{
-
- /*
- * This function performs the LU decomposition step. The coefficient
- * matrix 'amat' will be overwritten with the LU matrices in the
- * process.
- *
- * Parameters :
- * amat - This is the coefficient matrix, where the coefficients
- * are doubles. This function performs the row-major mapping
- * itself, declaring 'amat' to be a pointer rather than an array.
- * The calling routine can define the actual parameter to be
- * either.
- *
- * n - This parameter indicates the size of the current system
- * (nxn).
- *
- * numcols - This parameter indicates the actual defined column
- * dimension of the coefficient matrix. A bit reminiscent of
- * FORTRAN.
- */
-
- int i, j, k, n_minus_1 = n - 1;
- double dtmp1, *dptr1, *dptr2;
-
- /* To keep track of determinant sign due to row swaps*/
- int sign = 1;
-
- if (pivot != NULL) free(pivot);
- if ((pivot = (int *) malloc(n * sizeof(int))) == NULL) {
- fprintf(stderr, "Error in LUDecompose - malloc \n");
- return -2;
- }
-
- /* Initialize pivot vector*/
- for (i = 0; i < n_minus_1; i++) {
- *(pivot + i) = i;
- }
- *(pivot + n_minus_1) = n_minus_1;
-
- for (i = 0; i < n_minus_1; i++) {
- k = i;
- dptr1 = amat + i * numcols + i;
- for (j = i + 1; j < n; j++) {
- dptr2 = amat + j * numcols + i;
- if (fabs(*dptr2) > fabs(*dptr1) ) {
- dptr1 = dptr2;
- k = j;
- }
- }
-
- if (fabs(*dptr1) < SMALLEST_PIVOT) {
- fprintf(stderr, "Error in LUDecompose - matrix is singular\n");
- return -1;
- }
-
- /* Swap rows and update sign */
- if (k != i) {
- int temp = *(pivot + i);
- *(pivot + i) = *(pivot + k);
- *(pivot + k) = temp;
-
- /* Flip the sign of the determinant for each swap */
- sign *= -1;
-
- for (j = 0; j < n; j++) {
- dtmp1 = *(amat + i * numcols + j);
- *(amat + i * numcols + j) = *(amat + k * numcols + j);
- *(amat + k * numcols + j) = dtmp1;
- }
- }
-
- for (j = i+1; j < n; j++) {
- dtmp1 = *(amat + j * numcols + i) / *dptr1;
- *(amat + j * numcols + i) = dtmp1;
- for (k = i+1; k < n; k++) {
- *(amat + j * numcols + k) -= dtmp1 * *(amat + i * numcols + k);
- }
- }
- }
- *(pivot + n - 1) = sign;
- return 1;
-}
-
-void LUSolve(amat, b, n, numcols)
- double *amat, *b;
- int n, numcols;
- /*
- * This function solves a system of equations given a particular right
- * hand side vector 'b'. The coefficient matrix 'amat' must be
- * decomposed by the function LUDecompose prior to calling LUSolve.
- *
- * Parameters:
- * amat - This is the LU decomposition of the coefficient matrix.
- * The original coefficient matrix must first be passed to
- * the function LUDecompose before calling this routine.
- *
- * b - This parameter is the right hand side vector b. This
- * routine will compute the solution vector and return this
- * result in the parameter 'b'. NOTE - b WILL BE OVERWRITTEN
- *
- * n - This parameter indicates the size of the current system
- * (nxn).
- *
- * numcols - This parameter indicates the actual defined column
- * dimension of the coefficient matrix. A bit
- * reminiscent of FORTRAN.
- */
-{
- double dtmp1, *dptr1;
- int i, j, n_minus_1 = n - 1;
-
- /* Perform the row interchanges recorded in 'pivot' to vector b. */
- for (i = 0; i < n; i++) {
- /* Interchange element i and j. */
- j = *(pivot + i);
- if (j != i) {
- dtmp1 = *(b + i);
- *(b + i) = *(b + j);
- *(b + j) = dtmp1;
- }
- }
-
- /* Solve the system Ld = Pb for d. Vector d will overwritte b. */
- for (i = 0; i < n; i++) {
- dtmp1 = *(b + i);
- for (j = 0, dptr1 = amat + i*numcols; j < i; j++, dptr1++)
- dtmp1 -= *dptr1 * *(b + j);
- *(b + i) = dtmp1;
- }
-
- /* Solve the system Ux = d for x. Vector x will overwritte b (and d). */
- for (i = n_minus_1; i >= 0; i--) {
- dtmp1 = *(b + i);
- dptr1 = amat + i * numcols + n_minus_1;
-
- for (j = n_minus_1; j > i; j--) {
- dtmp1 -= *dptr1 * *(b + j);
- dptr1--;
- }
- *(b + i) = dtmp1 / *dptr1;
- }
-}
-
-double LUDeterminant(double *amat, int n, int numcols)
- /*
- * This function will compute the determinant of the matrix 'amat'.
- * 'amat' must have already been decomposed to LU form (use
- * LUDecompose).
- */
-{
- double det = 1.0;
- for (int i = 0; i < n; i++) {
- det *= *(amat + i * (numcols + 1));
- }
- return det * *(pivot + n - 1);
-}
+/* Function prototypes */
+static double factorial(int n);
+double Sin(double x);
+double Cos(double x);
+void matmul(double A[3][3], double B[3][3], double result[3][3]);
+void dot_product(double result[3], double A[3][3], double B[3]);
+void transpose(double A[3][3], double Ai[3][3]);
+void dispmat(double *A, int rows, int cols, bool transpose);
+int LUDecompose(double *amat, int n, int numcols);
+void LUSolve(double *amat, double *b, int n, int numcols);
+double LUDeterminant(double *amat, int n, int numcols);
+
+#endif /* UTILS_H */
