C语言代码实现PREBFORCE函数

#include <stdio.h> #include <math.h> #define NEH 100 #define L0H 100 #define NK 4

void PREBFORCE(double FH[NEH], double XH[L0H], double YH[L0H], double ZH[L0H], int IIIIH[NEH][NK], double BHF, double RAISE, int NEH, int L0H, int NK, int NCHE) { double XX[4], YY[4], ZZ[4], IX[4]; double DX21, DY21, DZ21, DX31, DY31, DZ31, DX41, DY41, DZ41; double XD1, XD2, XD3, AREAE1, AREAE2, AREAT = 0.0, BHF1, AREAE; int I, II, JJ; for(I = 0; I < NEH; I++) { FH[I] = 0.0; } for(II = 1; II <= NEH; II++) { for(JJ = 1; JJ <= NK; JJ++) { IX[JJ-1] = IIIIH[II-1][JJ-1]; XX[JJ-1] = XH[(int)(IX[JJ-1]-1)]; YY[JJ-1] = YH[(int)(IX[JJ-1]-1)]; ZZ[JJ-1] = ZH[(int)(IX[JJ-1]-1)]; } DX21 = XX[1] - XX[0]; DY21 = YY[1] - YY[0]; DZ21 = ZZ[1] - ZZ[0]; DX31 = XX[2] - XX[0]; DY31 = YY[2] - YY[0]; DZ31 = ZZ[2] - ZZ[0]; DX41 = XX[3] - XX[0]; DY41 = YY[3] - YY[0]; DZ41 = ZZ[3] - ZZ[0]; XD1 = DY21 * DZ31 - DY31 * DZ21; XD2 = DZ21 * DX31 - DZ31 * DX21; XD3 = DX21 * DY31 - DX31 * DY21; AREAE1 = 0.5 * sqrt(XD1 * XD1 + XD2 * XD2 + XD3 * XD3 + 1.0e-80); XD1 = DY41 * DZ31 - DY31 * DZ41; XD2 = DZ41 * DX31 - DZ31 * DX41; XD3 = DX41 * DY31 - DX31 * DY41; AREAE2 = 0.5 * sqrt(XD1 * XD1 + XD2 * XD2 + XD3 * XD3 + 1.0e-80); AREAT = AREAT + AREAE1 + AREAE2; } BHF1 = BHF - RAISE; if(BHF1 <= 0.0) { if(NCHE == 2) { printf('ERROR: Blank Holder force is unreasonable. '); fprintf(stderr, 'ERROR: Blank Holder force is unreasonable. '); } BHF1 = 0.0; } for(II = 1; II <= NEH; II++) { for(JJ = 1; JJ <= NK; JJ++) { IX[JJ-1] = IIIIH[II-1][JJ-1]; XX[JJ-1] = XH[(int)(IX[JJ-1]-1)]; YY[JJ-1] = YH[(int)(IX[JJ-1]-1)]; ZZ[JJ-1] = ZH[(int)(IX[JJ-1]-1)]; } DX21 = XX[1] - XX[0]; DY21 = YY[1] - YY[0]; DZ21 = ZZ[1] - ZZ[0]; DX31 = XX[2] - XX[0]; DY31 = YY[2] - YY[0]; DZ31 = ZZ[2] - ZZ[0]; DX41 = XX[3] - XX[0]; DY41 = YY[3] - YY[0]; DZ41 = ZZ[3] - ZZ[0]; XD1 = DY21 * DZ31 - DY31 * DZ21; XD2 = DZ21 * DX31 - DZ31 * DX21; XD3 = DX21 * DY31 - DX31 * DY21; AREAE1 = 0.5 * sqrt(XD1 * XD1 + XD2 * XD2 + XD3 * XD3 + 1.0e-80); XD1 = DY41 * DZ31 - DY31 * DZ41; XD2 = DZ41 * DX31 - DZ31 * DX41; XD3 = DX41 * DY31 - DX31 * DY41; AREAE2 = 0.5 * sqrt(XD1 * XD1 + XD2 * XD2 + XD3 * XD3 + 1.0e-80); AREAE = AREAE1 + AREAE2; FH[II-1] = BHF1 * AREAE / AREAT; } }