C INCLUDE 'CALPGM.INC' C miscellaneous subroutines: ULIB.FOR C C Copyright (C) 1989, California Institute of Technology C All rights reserved. U. S. Government Sponsorship under C NASA Contract NAS7-918 is acknowledged. C C Herbert M. Pickett, 20 March 1989 C SUBROUTINE CJJINI C set up log factorials REAL*8 FAC(0:1023),VN,V INTEGER NFAC COMMON /LNFACT/FAC SAVE /LNFACT/ C V=0 FAC(0)=0 FAC(1)=0 DO 500 NFAC=2,1023 VN=NFAC V=V+LOG(VN) 500 FAC(NFAC)=V RETURN END FUNCTION C3JJ(J1,J2,J3,M1,M2,M3) REAL*8 C3JJ INTEGER J1,J2,J3,M1,M2,M3 C C 3J COEFFICIENT CALCULATION C INTEGER JSUM,JM1,JM2,JM3,LA,LB,LC,LD,LE,IS,IFIN,I REAL*8 TRM,DELJJ,FAC(0:1023) COMMON /LNFACT/FAC SAVE /LNFACT/ C C3JJ=0 IF(M1+M2+M3.NE.0) RETURN TRM=DELJJ(J1,J2,J3,JSUM) JSUM=JSUM-1 IF(JSUM.LT.0)RETURN C CHECK FOR M1=M2=M3=0 CASE IF(M1.EQ.0) THEN IF(M3.EQ.0) THEN IS=JSUM/2 IF(JSUM.NE.IS+IS) RETURN LA=(JSUM-J1)/2 LB=(JSUM-J2)/2 LC=(JSUM-J3)/2 C3JJ=EXP(0.5*TRM+FAC(IS)-FAC(LA)-FAC(LB)-FAC(LC)) IF(IAND(IS,1).NE.0) C3JJ=-C3JJ RETURN ENDIF ENDIF I=J1+M1 IF(I.LT.0.OR.IAND(I,1).NE.0) RETURN JM1=I/2 TRM=TRM+FAC(JM1) JM1=J1-JM1 IF(JM1.LT.0) RETURN TRM=TRM+FAC(JM1) I=J2-M2 IF(I.LT.0.OR.IAND(I,1).NE.0) RETURN JM2=I/2 TRM=TRM+FAC(JM2) JM2=J2-JM2 IF(JM2.LT.0) RETURN TRM=TRM+FAC(JM2) I=J3-M3 IF(I.LT.0) RETURN JM3=I/2 TRM=TRM+FAC(JM3) JM3=J3-JM3 IF(JM3.LT.0) RETURN TRM=TRM+FAC(JM3) LA=JM1+J2-JSUM LB=JM2+J1-JSUM I=LA IF(I.LT.LB) I=LB IF(I.LT.0) I=0 LC=JSUM-J3 LD=JM1 LE=JM2 C3JJ=EXP(0.5*TRM-FAC(I)-FAC(I-LA)-FAC(I-LB) & -FAC(LC-I)-FAC(LD-I)-FAC(LE-I) ) IS=JSUM-J2+JM3+I IF(IAND(IS,1).NE.0) C3JJ=-C3JJ IFIN=LC IF(IFIN.GT.LD) IFIN=LD IF(IFIN.GT.LE) IFIN=LE IF(I.EQ.IFIN) RETURN LC=LC+1 LD=LD+1 LE=LE+1 IS=I+1 TRM=C3JJ DO 10 I=IS,IFIN TRM=-TRM*DBLE(LC-I)*DBLE(LD-I)*DBLE(LE-I)/ & (DBLE(I)*DBLE(I-LA)*DBLE(I-LB)) C3JJ=C3JJ+TRM 10 CONTINUE RETURN END FUNCTION C6JJ(J1,J2,J3,JJ1,JJ2,JJ3) REAL*8 C6JJ INTEGER J1,J2,J3,JJ1,JJ2,JJ3 C C 6J COEFFICIENT CALCULATION USING INTEGERS = 2 J C INTEGER L1,L2,L3,L4,LA,LB,LC,IS,IFIN,I REAL*8 TRM,DELJJ,FAC(0:1023) COMMON /LNFACT/FAC SAVE /LNFACT/ C C6JJ=0 TRM=DELJJ(J1,J2,J3,L1) IF(L1.EQ.0) RETURN I=L1 TRM=TRM+DELJJ(J1,JJ2,JJ3,L2) IF(L2.EQ.0) RETURN IF(I.LT.L2) I=L2 TRM=TRM+DELJJ(JJ1,J2,JJ3,L3) IF(L3.EQ.0) RETURN IF(I.LT.L3) I=L3 TRM=TRM+DELJJ(JJ1,JJ2,J3,L4) IF(L4.EQ.0) RETURN IF(I.LT.L4) I=L4 LA=L1+L4-J3 LB=L1+L3-J2 LC=L1+L2-J1 IS=I+1 C6JJ=EXP(0.5*TRM+FAC(I)-FAC(I-L1)-FAC(I-L2)-FAC(I-L3)-FAC(I-L4) & -FAC(LA-IS)-FAC(LB-IS)-FAC(LC-IS) ) IF(IAND(I,1).EQ.0) C6JJ=-C6JJ IFIN=LA IF(IFIN.GT.LB) IFIN=LB IF(IFIN.GT.LC) IFIN=LC IF(IS.EQ.IFIN) RETURN IFIN=IFIN-1 TRM=C6JJ DO 10 I=IS,IFIN TRM=-TRM*DBLE(I)*DBLE(LA-I)*DBLE(LB-I)*DBLE(LC-I)/ & ( DBLE(I-L1)*DBLE(I-L2)*DBLE(I-L3)*DBLE(I-L4) ) C6JJ=C6JJ+TRM 10 CONTINUE RETURN END FUNCTION DELJJ(JXA,JXB,JXC,JRET) REAL*8 DELJJ INTEGER JXA,JXB,JXC,JRET C C JXA,JXB,JXC=J+J C C ON RETURN JRET=J1+J2+J3+1 , UNLESS TRIANGLE RELATIONS NOT OBEYED C (IN WHICH CASE JRET=0) C C JRET IS ASSUMED TO BE LESS THAN 1024 C REAL*8 FAC(0:1023) INTEGER IA,IB,IC,IZ COMMON /LNFACT/FAC SAVE /LNFACT/ C JRET=0 DELJJ=0 IC=JXA+JXB+JXC IZ=IC/2 IF(IZ+IZ.NE.IC) RETURN IA=IZ-JXA IF(IA.LT.0) RETURN IB=IZ-JXB IF(IB.LT.0) RETURN IC=IZ-JXC IF(IC.LT.0) RETURN IZ=IZ+1 DELJJ=FAC(IA)+FAC(IB)+FAC(IC)-FAC(IZ) JRET=IZ RETURN END FUNCTION C9JJ(J11,J21,J31,J12,J22,J32,J13,J23,J33) REAL*8 C9JJ INTEGER J11,J21,J31,J12,J22,J32,J13,J23,J33 C C 9J COEFFICIENT CALCULATION C INTEGER JJ(0:8),LPERM,JSUM,JMIN,IXMIN,IX,J,JJTMP, & K1,K2,KK,JX REAL*8 C6JJ,VSQ,ONE PARAMETER (ONE=1.) C C9JJ=0 JJ(0)=J11 JJ(1)=J12 JJ(2)=J13 JJ(3)=J21 JJ(4)=J22 JJ(5)=J23 JJ(6)=J31 JJ(7)=J32 JJ(8)=J33 LPERM=0 C C FIND SMALLEST VALUE OF J C JSUM=JJ(0) JMIN=JSUM IXMIN=0 DO 106 IX=1,8 J=JJ(IX) JSUM=JSUM+J IF(J.LT.JMIN) THEN JMIN=J IXMIN=IX ENDIF 106 CONTINUE JSUM=JSUM/2 IF(IXMIN.GE.3) THEN C EXCHANGE COLUMN SO JMIN IS IN FIRST COLUMN IF(IXMIN.GE.6) THEN JX=6 ELSE JX=3 ENDIF IXMIN=IXMIN-JX DO 110 IX=0,2 JJTMP=JJ(IX) JJ(IX)=JJ(JX) JJ(JX)=JJTMP 110 JX=JX+1 LPERM=JSUM ENDIF IF(IXMIN.NE.0) THEN C EXCHANGE ROW SO JMIN IS IN FIRST ROW AND COLUMN IX=0 DO 120 J=1,3 JJTMP=JJ(IX) JJ(IX)=JJ(IX+IXMIN) JJ(IX+IXMIN)=JJTMP 120 IX=IX+3 LPERM=LPERM+JSUM ENDIF IF(JMIN.EQ.0) THEN C SPECIAL FOR J=0 IF(JJ(3).NE.JJ(6)) RETURN IF(JJ(1).NE.JJ(2)) RETURN VSQ=(JJ(3) + ONE)*(JJ(1) + ONE) C9JJ=C6JJ(JJ(8),JJ(7),JJ(3),JJ(4),JJ(5),JJ(1)) / SQRT(VSQ) LPERM=LPERM+ (JJ(1)+JJ(3)+JJ(5)+JJ(7))/2 ELSE C C FIND SMALLEST VALUE OF J NOT IN FIRST ROW OR COLUMN C JMIN=JJ(8) IXMIN=8 IF(JJ(4).LT.JMIN) THEN JMIN=JJ(4) IXMIN=4 ENDIF IF(JJ(5).LT.JMIN) THEN JMIN=JJ(5) IXMIN=5 ENDIF IF(JJ(7).LT.JMIN) THEN JMIN=JJ(7) IXMIN=7 ENDIF IF(IXMIN.LT.6) THEN C EXCHANGE COLUMN SO JMIN IS IN LAST COLUMN IXMIN=IXMIN+3 DO 130 IX=3,5 JJTMP=JJ(IX) JJ(IX)=JJ(IX+3) 130 JJ(IX+3)=JJTMP LPERM=LPERM+JSUM ENDIF IF(IXMIN.LT.8) THEN C EXCHANGE ROW SO JMIN IS IN LAST ROW AND COLUMN IX=1 DO 140 J=1,3 JJTMP=JJ(IX) JJ(IX)=JJ(IX+1) JJ(IX+1)=JJTMP 140 IX=IX+3 LPERM=LPERM+JSUM ENDIF K1= JJ(8)-JJ(0) KK= JJ(1)-JJ(5) IF(KK.LT.0) KK=-KK IF(K1.LT.KK) K1=KK KK= JJ(3)-JJ(7) IF(KK.LT.0) KK=-KK IF(K1.LT.KK) K1=KK K2= JJ(8)+JJ(0) KK= JJ(1)+JJ(5) IF(K2.GT.KK) K2=KK KK= JJ(3)+JJ(7) IF(K2.GT.KK) K2=KK IF(K1.GT.K2) RETURN DO 11 KK=K1,K2,2 11 C9JJ=C9JJ + DBLE(KK+1) & * C6JJ(JJ(0),JJ(1),JJ(2),JJ(5),JJ(8),KK) & * C6JJ(JJ(3),JJ(4),JJ(5),JJ(1),KK,JJ(7)) & * C6JJ(JJ(6),JJ(7),JJ(8),KK,JJ(0),JJ(3)) LPERM=LPERM+K1 ENDIF IF(IAND(LPERM,1).NE.0) C9JJ=-C9JJ RETURN END