PROGRAM PORCUPINE COMMON /BLOCK/ A(200),K REAL X(500), B(200), Y(200), XN(500) REAL T(200,200), S(500) REAL SUM, RUNIT, DY, DX,dis, TEGRAL, REZ, HN REAL,EXTERNAL::Y_OF_X INTEGER PGBEG, PGBAND, Z INTEGER JUNK, MODE CHARACTER*1 CH REAL PX, PY, PINAR, TUM, DEGER, SU !=======OPENING GRAPHICS DEVICE===================== IF (PGBEG(0,'?',1,1) .NE. 1) STOP !==========SETTING THE WINDOW PROPERTIES============= CALL PGPAGE CALL PGSVP(0.0,1.0,0.0,1.0) CALL PGSWIN(0.0,1.0,0.0,1.0) CALL PGBOX('bcts',0.1,5,'bcts',0.1,5) !=======COLLECTION OF DATA BY CURSOR================ PX = 0.5 PY = 0.5 MODE = 0 N=0 10 CONTINUE JUNK = PGBAND(MODE, 1, PX, PY, PX,PY,CH) WRITE (*, '(2F8.3,I4)') PX,PY,ICHAR(CH) IF (CH.EQ.'/') GOTO 20 CALL PGPT1(PX, PY,4) N = N + 1 X(N) = PX Y(N) = PY GOTO 10 !=======OBTAINING THE MAX. AND MIN. VALUES OF X=== 20 CONTINUE !=======The Maximum=============================== KB = 1 DO IF (KB == N ) EXIT DO I = 1,KB IF ( X(I) > X(I+1) ) THEN SU = X(I) X(I) = X(I+1) X(I+1) = SU ENDIF ENDDO KB = KB + 1 ENDDO AMAX = X(N) !=======The Minimum=============================== KA = N - 1 DO IF (KA == 0) EXIT DO I = 1,KA IF ( X(I) > X(I+1) ) THEN SU = X(I) X(I) = X(I+1) X(I+1) = SU ENDIF ENDDO KA = KA - 1 ENDDO AMIN = X(1) !=======OBTAINING THE COEFFICIENT MATRIX========== PRINT *,'ENTER THE ORDER' READ *,K DO I = 1,K DO J = 1,K SUM = 0 DO L = 1,N SUM = SUM + (X(L)**(J+I-2)) ENDDO T(I,J) = SUM ENDDO ENDDO !=======OBTAINING THE R.H.S MATRIX================= DO I = 1 , K SUM = 0 DO L = 1 , N SUM = SUM + Y(L)*(X(L)**(I-1)) ENDDO B(I) = SUM ENDDO !==========CALCULATION OF A1,A2,...,AK============== CALL GAUSS (T, B, A, K) !=======DRAWING THE K.th ORDER APPROXIMATION======== CALL PGSCI (2) 100 CALL PGFUNX (Y_OF_X , 1000 , AMIN, AMAX, 1) CALL PGSCI (7) PRINT *,'ENTER THE PORCUPINE NUMBER' READ *, Z PRINT *,'ENTER RUNIT' READ *, RUNIT !=======( K = ORDER )============================== !=======( Z = PORCUPINE NUMBER )=================== !=======ARRANGING POINTS TO PLACE PORCUPINES ON==== TUM = TRAPEZ (AMAX) EKSIK = TRAPEZ (AMIN) TUM = TUM - EKSIK dis = TUM / (Z+1) DO IC = 1,Z S(IC) = IC * dis ENDDO PRINT *,'Z 1. DEGER :',Z DO IB = 1 , Z XN(IB) = AMIN 300 XN ( IB ) = XN ( IB) + 0.001 AL = XN ( IB ) DEGER = TRAPEZ ( AL ) - TRAPEZ( AMIN ) REZ = S( IB ) - DEGER PINAR = ( REZ / S(IB) ) * 100 IF ( PINAR > 1.0 ) THEN GOTO 300 ENDIF IF ( REZ < 0 ) THEN XN ( IB )= AL ENDIF ENDDO PRINT *,'Z 2. DEGER :',Z !=======DRAWING PORCUPINES======================== DO L = 1,Z DY = ABS( ( ( Y_DER_TWO ( XN(L) ) * (RUNIT**2) )/ : ( 1.0 + Y_DER_ONE (XN(L) )**2)**2) ) DX = ABS( Y_DER_ONE ( XN(L) ) * DY) IF ( Y_DER_TWO ( XN(L) ) > 0 ) THEN YENX = XN ( L ) + SIGN ( DX , Y_DER_ONE ( XN ( L ) ) ) YENY = Y_OF_X ( XN ( L ) ) - DY ELSE IF ( Y_DER_TWO ( XN ( L ) ) < 0 ) THEN YENX = XN(L) - SIGN ( DX, Y_DER_ONE ( XN(L) ) ) YENY = Y_OF_X ( XN(L) ) + DY ELSE ENDIF ENDIF OBEZ = Y_OF_X ( XN ( L ) ) CALL PGMOVE ( XN ( L ) , OBEZ ) CALL PGDRAW ( YENX , YENY ) PRINT *,'Z L. DEGER :',L ENDDO PRINT *,'Z 3. DEGER :',Z PRINT *,'====WANT A NEW PORCUPINE CONFIGURATION?====' PRINT *,' ----PRESS 1 FOR YES ,2 FOR NO ,PLEASE---- ' READ *,NUM IF (NUM == 1)THEN CALL PGERAS GOTO 100 ELSE IF (NUM == 2) THEN GOTO 200 ENDIF ENDIF 200 CALL PGEND CONTAINS !=======FUNCTION FOR THE FIRST DERIVATIVE OF THE K.th ORDER POLYNOMIAL===== FUNCTION Y_DER_ONE (P) REAL Y_DER_ONE,P,SUM SUM = 0 DO I = 2,K SUM = SUM + (I-1) * A(I) * (P**(I-2)) ENDDO Y_DER_ONE = SUM END FUNCTION Y_DER_ONE FUNCTION FACT(X) REAL FACT INTEGER X,Z FACT = 1 DO Z = 2,X FACT = FACT * Z ENDDO END FUNCTION FACT !=======OBTAINING THE LENGTH OF THE CURVE================== FUNCTION TRAPEZ(B) REAL A,B,HN,XNO,TEGRAL,TRAPEZ A = 0.0 HN = (B-A)/250 TEGRAL = 0 DO IL = 1,249 XNO = A + ( IL * HN ) TEGRAL = TEGRAL+FONK(XNO) ENDDO TEGRAL = 0.5 *( FONK (A) + FONK (B) ) + TEGRAL TRAPEZ = TEGRAL * HN END FUNCTION TRAPEZ FUNCTION FONK(X) REAL FONK REAL,INTENT(IN)::X SUM = 0 DO I = 2 , K SUM = SUM + (I-1) * A(I) * (X**(I-2)) ENDDO FONK = SQRT ( (SUM**2) + 1 ) END FUNCTION FONK !=======FUNCTION FOR THE SECOND DERIVATIVE OF THE K.th ORDER POLYNOMIAL==== FUNCTION Y_DER_TWO ( P ) REAL Y_DER_TWO, P, SUM SUM = 0 DO I = 3,K SUM = SUM + FACT(I-1) * A(I) * (P**(I-3)) ENDDO Y_DER_TWO = SUM END FUNCTION Y_DER_TWO !=======GAUSS SUBROUTINE========================= SUBROUTINE GAUSS (a,b,r,n) REAL a(200,200),f(200,200),b(200),r(200) REAL s INTEGER x,z DO i = 1,n f(i,i) = 1.0 DO j = i+1,n f(i,j) = 0.0 ENDDO ENDDO DO 50 z = 1 , n - 1 DO x = z + 1 , n f(x,z) = a(x,z)/a(z,z) ENDDO DO 50 i = z+1 , n DO j = z,n a(i,j) = a(i,j)-(f(i,z)*a(z,j)) ENDDO b(i) = b(i)-f(i,z)*b(z) 50 CONTINUE r(n) = b(n)/a(n,n) DO i = n-1,1,-1 s = 0.0 DO j = i+1,n s = s + a(i,j) * r(j) ENDDO r(i) = 1/a(i,i) * (b(i)-s) ENDDO END SUBROUTINE GAUSS END PROGRAM PORCUPINE FUNCTION Y_OF_X ( X ) COMMON /BLOCK/ A(200),K REAL Y_OF_X REAL,INTENT(IN)::X SUM = 0 DO I = 1,K SUM = SUM + ( A(I) *( X**(I-1) ) ) ENDDO Y_OF_X = SUM END FUNCTION Y_OF_X