      PROGRAM WOOF
C Program to plot the frequency response of a loudspeaker enclosure
C given the speaker's parameters. From IEEE Trans. Electroacoustics.
      DIMENSION RESP(49),ZEL(49),JY(4),JZ(4),J(6),MATRIX(49,29),MK(3)
     1,J1(6)
      REAL M,M1,M2,M2OPT
      DATA A1,V,FAR,S1,BL,R,Q0,M2/.051,.45,40.6,1777.,7.9,5.5,
     15.09,.013/
      DATA MK(1),MK(2),MK(3)/1H ,1H+,1H./
      DATA JY/-10,-20,-30,-40/
      DATA JZ/30,20,10,0/
      DATA J/0,40,80,120,160,200/
      DATA J1/0,10,20,30,40,50/
    6 WRITE (3, 90)
   90 FORMAT('1LOUDSPEAKER ENCLOSURE ANALYSIS PROGRAM.'/
     1'0For information about this program contact'/
     2' Trevor MARSHALL, SYSOP, T.O. Tech RBBS'/
     3'           805-492-5472, (voice)805-492-3693')
       WRITE (3, 567)
  567 FORMAT(' '/' For HELP type 1,otherwise 0'/)
      READ (3,568) IHELP
  568 FORMAT(I1)
      IF(IHELP.EQ.1)GO TO 676
      WRITE (3, 452)
  452 FORMAT(/
     1'0WHAT IS THE EFFECTIVE PISTON AREA IN SQUARE METRES?')
      READ (3,1)A1
      WRITE (3, 91)
   91 FORMAT(' ADIABATIC ENCLOSURE VOLUME IN CUBIC METRES?'/
     1'(SET TO 99999 IF AN INFINITE BAFFLE IS IN USE)')
      READ (3,1)V
      WRITE (3, 92)
   92 FORMAT(' WOOFER FREE AIR RESONANCE?')
      READ (3,1)FAR
      WRITE (3, 93)
   93 FORMAT(' SUSPENSION STIFFNESS IN NEWTON/METRE?')
      READ (3,1)S1
      WRITE (3, 94)
   94 FORMAT(' B L PRODUCT IN WEBERS/METER?')
      READ (3,1)BL
      WRITE (3, 95)
   95 FORMAT(' VOICE COIL RESISTANCE?')
      READ (3,1)R
      WRITE (3, 96)
   96 FORMAT(' SUSPENSION Q FACTOR?')
      READ (3,1)Q0
      WRITE (3, 97)
   97 FORMAT(' VENT AIR MASS IN KILOGRAMS'/
     1' (SET TO 9999 IF ENCLOSURE IS SEALED)')
      READ (3,1)M2
    1 FORMAT(F10.0)
      IF(BL.EQ.0.)STOP
  676 CONTINUE
C      WRITE (3,  999)
 999  FORMAT(1H1)
      WRITE (3,  10)
  10  FORMAT(21X,26HWOOFER PERFORMANCE PLOTTER  /20X,28(1H'),/'0THIS PRO
     1GRAM CALCULATES IN TWO FREQUENCY RANGES,'/' 0 TO 200 HZ IN 4 HZ STE
     2PS '/' AND 0 TO 50 HZ IN 1 HZ STEPS.')
      WRITE (3, 20)A1,V,FAR,S1,BL,R,Q0,M2
   20 FORMAT(4X,24HSYSTEM INPUT DATA ARE...  /
     228H EFFECTIVE PISTON AREA       ,F6.3,17H SQUARE METRES   /
     328H ADIABATIC ENCLOSURE VOLUME  ,F6.3,17H CUBIC METRES    /
     428H WOOFER FREE AIR RESONANCE   ,F6.1,17H HERTZ           /
     528H SUSPENSION STIFFNESS        ,F6.0,17H NEWTON/METRE    /
     628H B L PRODUCT                 ,F6.1,17H WEBER/METRE     /
     728H VOICE COIL RESISTANCE       ,F6.1,17H OHMS            /
     828H SUSPENSION Q FACTOR         ,F6.2,17H                 /
     928H REFLEX VENT AIR MASS        ,F6.3,' KILOGRAM')
      IF(M2.GT.999.) WRITE (3, 2222)
2222  FORMAT(' THIS ENCLOSURE WILL BE TREATED AS FULLY SEALED (INFINITELY
     1Y SMALL VENT).  ')
CONCLUDES INPUT CHECK
CALCULATE PARAMETERS
      M1=S1/(6.283*FAR)**2
      S2=143000.*A1**2/V
      M2OPT=M1*S2/S1
      REM=BL**2/R
      FH=0.159*SQRT(S2/M2)
      ASYMP=0.0552*REM*(A1/M1)**2
      X1=SQRT(S1*M1)
      R1=X1/Q0
      Q=X1/(R1+REM)
      S=S2/S1
      M=M2/M1
      A=Q**(-2.)-2.-2.*S-2.*S/M
      B=1.+2.*S+S*S+4.*S/M+2.*S*S/M+S*S/M/M-2.*S/Q/Q/M
      C=S*S/Q/Q/M/M-2.*S/M-2.*S*S/M-2.*S*S/M/M
      D=S*S/M/M
      WRITE (3,  30)M1,S2,REM,FH,M2OPT,ASYMP,Q,S,M,A,B,C,D
  30  FORMAT(4X,43HFROM THESE DATA, CONSTANT PARAMETERS ARE... /
     228H MASS OF CONE AND AIR LOAD    ,F6.3,17H KILOGRAM         /
     328H ENCLOSURE AIR STIFFNESS     ,F6.0,17H NEWTON/METRE      /
     428H ELECTRODYNAMIC DRAG         ,F6.2,17H NEWTON SEC/METRE  /
     528H HELMHOLTZ FREQUENCY         ,F6.1,17H HERTZ             /
     628H SUGGESTED VENT AIR MASS     ,F6.3,17H KILOGRAM          /
     728H ASYMPTOTIC EFFICIENCY     ,F6.3,17H PERCENT             /
     8' DIAGNOSTIC PARAMETERS Q,S,M,A,B,C,D,ARE'7(1X,F5.2))
COMMENCE EVALUATING RESPONSE AND IMPEDANCE
      DO 998 ISTEP=1,2
C      IF(ISTEP.EQ.2)WRITE (3, 999)
      DO100 I=1,49
      IF(ISTEP.EQ.1)FR=4*I
      IF(ISTEP.EQ.2)FR=I
      W=6.283*FR
      RRAD=.022*(FR*A1)**2
      R2=RRAD
	FR=W*M2
	FI=FR-S2/W
	CALL CDIV(R2,FR,R2,FI,FR1,FI1)
	ZR=(R1+RRAD)
	ZI=(W*M1-S1/W)
	CALL CMULT(0.,1.,FR1,FI1,FR2,FI2)
	ZR=ZR-FR2*S2/W
	ZI=ZI-(FI2*S2/W)
	ZR2=ZR+REM
	RESP(I)=RRAD*REM*CABS(FR1,FI1)**2/CABS(ZR2,ZI)**2
	CALL CDIV(REM,0.,ZR,ZI,ZR1,ZI1)
	ZR1=ZR1+1.
	ZEL(I)=R*CABS(ZR1,ZI1)
100	CONTINUE
CLEAR AND FILL DISPLAY MATRIX
      DO 200 IFREQ=1,49
      DO150 LEVEL=1,29
 150  MATRIX(IFREQ,LEVEL)=MK(1)
      IY=40.5+10.*ALOG10(RESP(IFREQ))
      IF((IY.LT.1).OR.(IY.GT.29))GO TO 162
  162 MATRIX(IFREQ,IY)=MK(2)
 170  IZ=0.5+ZEL(IFREQ)
      IF((IZ.LT.1).OR.(IZ.GT.29))GO TO 200
 172  MATRIX(IFREQ,IZ)=MK(3)
 200  CONTINUE
C   NOW WRITE  DISPLAY MATRIX
      WRITE (3,  600)
 600  FORMAT('1',8X,'ABSOLUTE RESPONSE(+)IN DB,AND IMPEDANCE(.)IN OHMS')

      WRITE (3, 601)
 601  FORMAT(6(9X,1H+))
      WRITE (3, 602)JY(1),JZ(1)
 602  FORMAT(I9,51(1H+),I3)
      DO 670 LOSS=1,29
      IY=30-LOSS
      IF(LOSS-10)615,616,615
 615  IF(LOSS-20) 617,618,617
 617  WRITE (3,  620)(MATRIX(IFREQ,IY),IFREQ=1,49)
 620  FORMAT(9X,1H+,49A1,1H+)
      GO TO 670
 616  WRITE (3,  630)JY(2),(MATRIX(IFREQ,IY),IFREQ=1,49),JZ(2)
 630  FORMAT(I9,1H+,49A1,1H+,I3)
      GO TO 670
 618  WRITE (3,  630)JY(3),(MATRIX(IFREQ,IY),IFREQ=1,49),JZ(3)
  670 CONTINUE
      WRITE (3, 602)JY(4),JZ(4)
      WRITE (3,  601)
      IF(ISTEP.EQ.2)WRITE (3,  680)J1
      IF(ISTEP.EQ.1)WRITE (3,  680)J
 680  FORMAT(6(7X,I3))
      WRITE (3, 690)
 690  FORMAT(25X,17HFREQUENCY,  HERTZ   ///)
 998  CONTINUE
      WRITE (3, 1000)
 1000 FORMAT(' TYPE 0 TO FINISH, 1 FOR ANOTHER RUN')
      READ (3,1001)IRUN
 1001 FORMAT(I1)
      IF(IRUN.EQ.1)GO TO 6
      STOP
	END
	FUNCTION CABS(X,Y)
	CABS=SQRT(X*X+Y*Y)
	RETURN
	END
	SUBROUTINE CMULT(OPR1R,OPR1I,OPR2R,OPR2I,RLTR,RLTI)
	RLTR=OPR1R*OPR2R-OPR1I*OPR2I
	RLTI=OPR1R*OPR2I+OPR1I*OPR2R
	RETURN
      END
	SUBROUTINE CDIV(OPR1R,OPR1I,OPR2R,OPR2I,RLTR,RLTI)
	DIV9=OPR2R*OPR2R+OPR2I*OPR2I
	RLTR=OPR1R*OPR2R-OPR1I*(0.-OPR2I)
	RLTI=OPR1R*(0.-OPR2I)+OPR1I*OPR2R
	RLTR=RLTR/DIV9
	RLTI=RLTI/DIV9
	RETURN
	END
