PROGRAM OBR5

C Modifikace vyrazu pro rudy posuv.
C
C Modeluje obraz akrecniho disku kolem cerne diry
C a bere v uvahu prvni a druhy rad

      PARAMETER( NFITS = 100 , NFITS2 = 10000 )

      INTEGER L,K,FOPEN,NOBR(NFITS,NFITS)
      REAL R,DFI,FI,FIN,BFI,INCL,P,SUCFAI,FE,B,X,Y,T,FIC
      REAL ES,LS,G1,G2,SINCL,COSFI
      REAL OBR(NFITS,NFITS)
      EXTERNAL SUCFAI,FE
      COMMON /DAT/ R,FIN                             
      SAVE /DAT/
                                   
      DATA INCL/89.0/,OBR/NFITS2*0.0/,NOBR/NFITS2*0/
                            
      INCL = INCL/57.29578                      
      SINCL = SIN(INCL)
      DO R = 6,40,0.01
C Bod odkud neni navratu
        T = SQRT(6.0/R + 1.0)
        FIC = LOG(0.26794919*(1.7320508 + T)/(1.7320508 - T))
        P = 5.0
        DO DFI = 0,179,0.05
          FI = DFI/57.29578
C paprsky prvniho radu
          COSFI= COS(FI)
          FIN = ACOS( SINCL*COSFI )
C          IF( FIN.LT.FIC ) GOTO 50
          BFI = ASIN(SIN(FI)/SQRT(1.0-(SINCL*COSFI)**2))
          IF( FI.GT.1.5707963 ) BFI = 3.1415927 - BFI
          IF( SUCFAI(FE,P,0.1,5E-3,50).GT.1E-2 )THEN
            PRINT *,'!>50',int(DFI),P
            GOTO 50
             P = R
          ENDIF
          B = SQRT( P**3/(P - 2.0) )
          X = B*SIN(BFI)
          Y = -B*COS(BFI)       
c          PRINT *,X,Y
c         b = R*SIN(FIN)
c          X = B*SIN(BFI)
c          Y = -B*COS(BFI)       
c          PRINT *,X,Y
c          PRINT  *,-X,Y
          LS = SQRT(R**2/(R - 3.0))
          ES = SQRT((R - 2.0)**2/R/(R - 3.0))
          G2 = ES/(1.0 - 2.0/R)
          T = B*LS*SINCL*SIN(BFI)/R**2
          G1 = SQRT(1.0 - 2.0/R)*(G2 + T)
          G2 = SQRT(1.0 - 2.0/R)*(G2 - T) 
          G1 = 1.0/G1**4
          G2 = 1.0/G2**4
          L = (X + 50.0)*(NFITS/100.0)
          K = (Y + 50.0)*(NFITS/100.0)
          IF( L.GE.1 .AND. L.LE.NFITS )THEN
            IF( K.GE.1 .AND. K.LE.NFITS )THEN
C              OBR(L,K) = 900.0
C              OBR(NFITS - L + 1,K) = 900.0
              OBR(L,K) = OBR(L,K) + G1
              OBR(NFITS-L+1,K) = OBR(NFITS-L+1,K) + G2
              NOBR(L,K) = NOBR(L,K) + 1
              NOBR(NFITS-L+1,K) = NOBR(NFITS-L+1,K) + 1
            ENDIF
          ENDIF
C paprsky druheho radu
          FIN = 6.2831853 - FIN  
          IF( SUCFAI(FE,P,0.1,5E-3,50).GT.1E-2 )THEN
            PRINT *,'!>50',int(DFI),P
            GOTO 50
             P = R
          ENDIF
          B = SQRT( P**3/(P - 2.0) )
          X = B*SIN(BFI)
          Y = B*COS(BFI)       
c          PRINT *,X,Y
c         b = R*SIN(FIN)
c          X = B*SIN(BFI)
c          Y = -B*COS(BFI)       
c          PRINT *,X,Y
c          PRINT  *,-X,Y
          G2 = ES/(1.0 - 2.0/R)
          T = B*LS*SINCL*SIN(BFI)/R**2
          G1 = G2 + T
          G2 = G2 - T 
          G1 = 1.0/G1**4
          G2 = 1.0/G2**4
          L = (X + 50.0)*(NFITS/100.0)
          K = (Y + 50.0)*(NFITS/100.0)
          IF( L.GE.1 .AND. L.LE.NFITS )THEN
            IF( K.GE.1 .AND. K.LE.NFITS )THEN
C              OBR(L,K) = 900.0
C              OBR(NFITS - L + 1,K) = 900.0
              OBR(L,K) = OBR(L,K) + G1
              OBR(NFITS-L+1,K) = OBR(NFITS-L+1,K) + G2
              NOBR(L,K) = NOBR(L,K) + 1
              NOBR(NFITS-L+1,K) = NOBR(NFITS-L+1,K) + 1
            ENDIF
          ENDIF
50      ENDDO
      ENDDO
          
      T = 0.0
      DO L = 1,NFITS
        DO K = 1,NFITS
          IF( NOBR(K,L).NE.0 ) OBR(K,L) = OBR(K,L)/NOBR(K,L)
          T = MAX(OBR(K,L),T)
        ENDDO
      ENDDO
      
      IF( FOPEN(2,'d.fts',NFITS).NE.0 ) STOP 002 
      
      DO L = 1,NFITS
        DO K = 1,NFITS      
          CALL PUTPIX(2,int(1E4*OBR(K,L)/T))
        ENDDO
      ENDDO
      
      CALL FCLOSE(2)
      
      END
      
      REAL FUNCTION FE(P)
      REAL P
      
      REAL R,FIN,Q,K2,CHI2,CHI2N,ALPHA,F,ELIPT1,K1,K,FIT
      REAL R0
      COMMON /DAT/ R,FIN
      SAVE /DAT/,FIT,R0
      DATA R0/0.0/
                
      IF( R0.NE.R )THEN     
        R0 = R
C Vypocet uvrate
        Q = SQRT( (R - 2.0)*(R + 6.0) )
        K2 = 0.5*(Q - R + 6.0)/Q
        CHI2N = ACOS( SQRT(2.0/Q/K2) )
        K = SQRT(K2)
        ALPHA = ASIN(SQRT(K))
        FIT = 2.0*SQRT(R/Q)*(K1(K) - ELIPT1(CHI2N,ALPHA))
      ENDIF
      IF( P.LE.3.0 .OR. P.GT.R )THEN
        FE = 1E10
        RETURN
      ENDIF  
      Q = SQRT( (P - 2.0)*(P + 6.0) )
      K2 = 0.5*( Q - P + 6.0 )/Q
      CHI2N = ACOS( SQRT(2.0/Q/K2) )
      CHI2  = ACOS( SQRT(2.0*(1.0 - P/R)/Q/K2) )
      K = SQRT(K2)
      ALPHA = ASIN(SQRT(K))
      IF( FIN.LT.FIT )THEN
        F = 2.0*SQRT(P/Q)*(ELIPT1(CHI2,ALPHA)-ELIPT1(CHI2N,ALPHA))
      ELSE
        F = 2.0*SQRT(P/Q)*
     *     (2.0*K1(K) - ELIPT1(CHI2,ALPHA) - ELIPT1(CHI2N,ALPHA))
      ENDIF
      FE = (FIN - F)**2
c      IF( ABS(FIN-0.885).LE.0.005.or.ABS(FIN-2.8955).LE.0.005 )THEN
c         PRINT '(10F7.3)',P,Q,SQRT(K2),CHI2*57.3,CHI2N*57.3,F*57.3,
c     *          FIN*57.3,FIT*57.3
c      ENDIF
      END