FUNCTION LOGPOT(X)
        REAL*4 X, LOGPOT
        REAL*4 YM, A, B
        
        DATA  YM /0.2/
        
        B = ((1. / YM) - 1.) ** 2.
        
        A = 1. / (B - 1.)

        LOGPOT = A * (B ** X) - A
        END FUNCTION


        FUNCTION LOGPRO(V, T1, T2)
        REAL*4 V, T1, T2, LOGPRO
        REAL*4 LOGPOT
        REAL*4 LV, LT1, LT2

        LV = LOGPOT(V)
        LT1 = LOGPOT(T1)
        LT2 = LOGPOT(T2)

        LOGPRO = (LV - LT1) / (LT2 - LT1)
        END FUNCTION

        FUNCTION PAR(A, B)
        COMPLEX*8 A, B, PAR

        PAR = 1. / ((1. / A) + (1. / B))
        END FUNCTION

        FUNCTION CAP(C, F)
        REAL*4 C, F, PI
        COMPLEX*8 CAP

        DATA PI /3.14159265358979/
        
        CAP = (0, -1) / ( C * 2 * PI * F)
        END FUNCTION

        FUNCTION GAIN(F, K)
        REAL*4 F, K, GAIN

        REAL*4 LOGPOT
        REAL*4 LOGPRO
        COMPLEX*8 PAR
        COMPLEX*8 CAP

        REAL*4 TE, TD, RAC
        REAL*4 RAE, RAD, RED, RDC, TOT
        COMPLEX*8 LD, LE

        COMPLEX*8 ERAE, ERAD, ERED, ERAC, ERDC, ETOT
        COMPLEX*8 PE, PD, PA, PED, PDA, G
        REAL*4 V

        DATA TE /0.4/
        DATA TD /0.6/
        DATA RAC /470.E3/

        LD = 18.E3 + CAP( 47.E-9, F)
        LE = 4.7E3 + CAP( 250.E-9, F)

        
        RAE = RAC * LOGPOT(TE)
        RAD = RAC * LOGPOT(TD)

        RED = RAD - RAE
        RDC = RAC - RAD

        TOT = RAE + RED + RDC

        WRITE(*,*) "TRACK RESISTANCES: ", RAE, "+", RED, "+", 
     C      RDC, "=", TOT

        ERAE = PAR(RAE * (1, 0.), LE)
        ERAD = PAR(ERAE + RED, LD)
        ERED = ERAD - ERAE
        ERAC = ERAD + RDC
        ERDC = ERAC - ERAD

        ETOT = ERAC + 150000. + CAP(15.E-9, F) + 
     C      PAR((1100., 0.), CAP(120.E-9, F))

        WRITE(*, *) "EFFECTIVE TOTAL IMPEDANCE", ERAE, "+",
     C      ERED, "+", ERDC, "=", ERAC

        PE  = ERAE / ETOT
        PD  = ERAD / ETOT
        PA  = ERAC / ETOT
        PED = PD - PE
        PDA = PA - PD
    
        V = K / 10.

        IF (V.LT.TE) THEN
          G = PE * LOGPRO(V, 0., TE)
        ELSEIF (V.LT.TD) THEN
          G = PE + (PED * LOGPRO(V, TE, TD))
        ELSE
          G = PD + ( PDA * LOGPRO( V, TD, 1. ) )
        ENDIF

        GAIN = ABS(G)
        END FUNCTION


        PROGRAM MAIN
        INTEGER*4 K
        REAL*4 F, G, GMAX
        CHARACTER*13 FN
        CHARACTER*1 FNA
        INTEGER*1 INA
        DIMENSION FNA(13), INA(13)
        EQUIVALENCE(FN,FNA,INA)
        DATA FNA /'V','O','L','U','M','E','-','0','0','.','D','A','T'/


        OPEN (UNIT = 10, FILE = 'VOLUME.PLT', ACCESS = 'SEQUENTIAL', 
     C      STATUS = 'UNKNOWN')

        WRITE(10, *) 'set term png size 1280,1024'
        WRITE(10, *) 'set output "volume.png"'
        WRITE(10, *) 'set logscale x'
        WRITE(10, *) 'set logscale y'
        WRITE(10, *) 'set yrange [:1.1]'
        WRITE(10, *) 'set xlabel "frequency/Hz"'
        WRITE(10, *) 'set ylabel "gain"'
        WRITE(10, *) 'plot \';

        DO 40 K = 1, 10
 
        FNA(8) = '0'
        INA(8) = INA(8) + (K / 10)
        FNA(9) = '0'
        INA(9) = INA(9) + MOD(K, 10)

        WRITE(10, *) '    "', FN, '" using 1:2 with lines title "',
     C      K,'", \'

        OPEN (UNIT = 11, FILE = FN, ACCESS = 'SEQUENTIAL', 
     C      STATUS = 'UNKNOWN')


        GMAX = 0.

        F = 20.
10      G = GAIN(F, 1. * K)
        IF (G.GT.GMAX) THEN
          GMAX = G
        END IF
        F = F * 1.1
        IF (F - 20.E3) 10,10,20
20      F = 20.
30      G = GAIN(F, 1. * K)
        WRITE(11, *) F, G, G / GMAX
        F = F * 1.1
        IF (F - 20.E3) 30,30,40
40      CONTINUE

        CLOSE(10)

        WRITE(11,*) ''
        CLOSE(11)

        CALL SYSTEM('gnuplot VOLUME.PLT')

        END PROGRAM