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