       PROGRAM POTEN1_3D
C++::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
C
C  This program calculates a Schwarzschild's model
C  of a triaxial galaxy using De Zeeuw's approximation (1983)
C   and projects the density on one of the principal planes
C
C JLP Version of 10-10-86
C--::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
C
	PARAMETER (IDIM=1000)
	PARAMETER (IDIM1=50)
	IMPLICIT DOUBLE PRECISION (A-H,O-Z)
	REAL*4 IMAGE(IDIM,IDIM)
	REAL*8 PHI(0:IDIM1,0:IDIM1,0:IDIM1)
	REAL*8 RHO(0:IDIM1,0:IDIM1,0:IDIM1)
	CHARACTER NAME*40,COMMENTS*80,ANS*1,AGAIN*1
	COMMON/POT1/PI,QQ,PP,PASX,PASY,PASZ
C
	PI=3.14159265358979323846
C
C GRAV=6.670 E-8 in dyn.cm2.g-2  (C.G.S.), but like QUINN
C we take GRAV=1.
	GRAV=1.D0
 
	PRINT *,' PROGRAM POTEN VERSION OF JUNE 8th 1986'
 
C Inquire the format of the files :
	CALL JLP_INQUIFMT
 
	OPEN(1,FILE='poten.log',STATUS='unknown',
     1	ACCESS='SEQUENTIAL')
C
C  Input of the parameters
C
10	FORMAT(A)
 
	PRINT *,' NUMBER OF POINTS FOR THE MODEL (LESS THAN 50)?'
	READ(5,*) NPTS
1000	WRITE (6,335)
335	FORMAT(' Enter P, Q : ',$)
	READ (5,*) PP,QQ
	PRINT *,' ENTER PASX, PASY, PASZ :'
	READ(5,*) PASX,PASY,PASZ
 
	WRITE(1,340),PP,QQ,PASX,PASY,PASZ
340	FORMAT(4X,'P =',F8.3,'  Q =',F8.3,5X,'PASX, PASY, PASZ :',
     1	3(F9.3,2X),/)
 
C Calculates C1 and C3
 
	PP2=PP*PP
	QQ2=QQ*QQ
	PP3=PP2*PP
	QQ3=QQ2*QQ
	ASYM2=2.D0*PP2*QQ2-PP2-QQ2
	SYM2=PP2*QQ2+PP2+QQ2
	SYM3=1.D0+PP3+QQ3
 
	D1=ASYM2/(14.D0*SYM2)
	D2=(1.D0-(3.D0/SYM3))/6.D0
	D3=(PP2-QQ2)/(28.D0*SYM2)
	D4=(PP3-QQ3)/(12.D0*SYM3)
 
C	C1=(4.D0/35.D0)*ASYM2/RRR
C	C3=(2.D0/35.D0)*(PP2-QQ2)/RRR
 
	C1=(D1**0.6)*(D2**0.4)
	C2=(D1/D2)**0.4
	C3=(D3**0.6)*(D4**0.4)
	C4=(D3/D4)**0.4
 
	CTE1=(4.D0/PI)/(4*PI*GRAV)
 
	WRITE(1,341) C1,C2,C3,C4
341	FORMAT(4X,'C1 =',E12.5,'  C2 =',E12.5,
     1	'C3 =',E12.5,'  C4 =',E12.5,/,/,11X,'X',
     1	13X,'Y',13X,'Z',10X,'PHI(X,Y,Z)',7X,'RHO(X,Y,Z)',/)
 
	DO 100 IX=0,NPTS
	X=FLOAT(IX)*PASX
	X2=X*X
 
	DO 100 IY=0,NPTS
	Y=FLOAT(IY)*PASY
	Y2=Y*Y
 
	DO 100 IZ=0,NPTS
	Z=FLOAT(IZ)*PASZ
	Z2=Z*Z
 
	RAD2=X2+Y2+Z2
	IF(RAD2.EQ.0.) RAD2=1.E-12
	RAD=DSQRT(RAD2)
 
C Calculating U, V and W
 
	URAD=(-1.D0/RAD)*DLOG(RAD+DSQRT(1+RAD2))+1.D0
	VRAD=(-1.D0*C1*RAD2)/((1.D0+C2*RAD2)**1.5)
	WRAD=(-1.D0*C3*RAD2)/((1.D0+C4*RAD2)**1.5)
	
C Remember :CTE1=(4.D0/PI)/(4*PI*GRAV)
C
	FRAD=1.D0/((1.D0+RAD2)**1.5)
	GRAD=3.D0*C1*C2*(7.D0+2.D0*C2*RAD2)*RAD2/
     1	((1.D0+C2*RAD2)**3.5)
	HRAD=3.D0*C3*C4*(7.D0+2.D0*C4*RAD2)*RAD2/
     1	((1.D0+C4*RAD2)**3.5)
 
 
	AA=(2.D0*Z2-X2-Y2)/(2.D0*RAD2)
	BB=3.D0*(X2-Y2)/RAD2
 
C
	PHI(IX,IY,IZ)=URAD-AA*VRAD+BB*WRAD
	RHO(IX,IY,IZ)=FRAD-AA*GRAD+BB*HRAD
 
	IF(IY.EQ.1.AND.IZ.EQ.1)THEN
	WRITE(1,342)X,Y,Z,PHI(IX,IY,IZ),RHO(IX,IY,IZ)
	ENDIF
342	FORMAT(2X,3(F12.3,2X),4X,E12.5,4X,E12.5)
 
100	CONTINUE
 
	NPL=2*NPTS+1
	NL=2*NPTS+1
 
	PRINT *,' DO YOU WANT TO STORE THE IMAGE OF THE DENSITY'
 
	PRINT *,' PROJECTED ON THE XY PLANE ?'	
	READ(5,10) ANS
	IF(ANS.EQ.'Y') THEN
	FACT=1.E+04
	CALL PROJECT(RHO,NPTS,IDIM1,IMAGE,IDIM,FACT,1)
        WRITE(6,*) 'Output file: '
        READ(5,10) NAME
	COMMENTS=' '
	CALL JLP_WRITEIMAG(IMAGE,NPL,NL,IDIM,NAME,COMMENTS)
	ENDIF
 
	PRINT *,' PROJECTED ON THE YZ PLANE ?'	
	READ(5,10) ANS
	IF(ANS.EQ.'Y') THEN
	FACT=1.E+04
	CALL PROJECT(RHO,NPTS,IDIM1,IMAGE,IDIM,FACT,2)
        WRITE(6,*) 'Output file: '
        READ(5,10) NAME
	COMMENTS=' '
	CALL JLP_WRITEIMAG(IMAGE,NPL,NL,IDIM,NAME,COMMENTS)
	ENDIF
 
	PRINT *,' PROJECTED ON THE ZX PLANE ?'	
	READ(5,10) ANS
	IF(ANS.EQ.'Y') THEN
	FACT=1.E+04
	CALL PROJECT(RHO,NPTS,IDIM1,IMAGE,IDIM,FACT,3)
        WRITE(6,*) 'Output file: '
        READ(5,10) NAME
	COMMENTS=' '
	CALL JLP_WRITEIMAG(IMAGE,NPL,NL,IDIM,NAME,COMMENTS)
	ENDIF
 
	PRINT *,' DO YOU WANT TO STORE THE IMAGE OF THE POTENTIAL?'
 
	PRINT *,' 1: IN THE PLANE X,Y ?'
	READ(5,10) ANS
	IF(ANS.EQ.'Y') THEN
	FACT=-1.
	CALL PLANE(PHI,NPTS,IDIM1,IMAGE,IDIM,FACT,1)
        WRITE(6,*) 'Output file: '
        READ(5,10) NAME
	COMMENTS=' '
	CALL JLP_WRITEIMAG(IMAGE,NPL,NL,IDIM,NAME,COMMENTS)
	ENDIF
 
	PRINT *,' 1: IN THE PLANE Y,Z ?'
	READ(5,10) ANS
	IF(ANS.EQ.'Y') THEN
	FACT=-1.
	CALL PLANE(PHI,NPTS,IDIM1,IMAGE,IDIM,FACT,2)
        WRITE(6,*) 'Output file: '
        READ(5,10) NAME
	COMMENTS=' '
	CALL JLP_WRITEIMAG(IMAGE,NPL,NL,IDIM,NAME,COMMENTS)
	ENDIF
 
	PRINT *,' 1: IN THE PLANE Z,X ?'
	READ(5,10) ANS
	IF(ANS.EQ.'Y') THEN
	FACT=-1.
	CALL PLANE(PHI,NPTS,IDIM1,IMAGE,IDIM,FACT,3)
        WRITE(6,*) 'Output file: '
        READ(5,10) NAME
	COMMENTS=' '
	CALL JLP_WRITEIMAG(IMAGE,NPL,NL,IDIM,NAME,COMMENTS)
	ENDIF
 
C
999	CLOSE(1)
	PRINT *,' Output in "poten.log".'
	STOP
	END
C
C
 
	SUBROUTINE PROJECT(RHO,NPTS,IDIM1,IMAGE,IDIM,FACT,IOPT)
C****************************************************
C	Subroutine PROJECT to project the density on one
C	of the planes XY (IOPT=1), YZ (IOPT=2) or ZX (IOPT=3)
C********************************************************
	IMPLICIT DOUBLE PRECISION (A-H,O-Z)
	REAL*8 RHO(0:IDIM1,0:IDIM1,0:IDIM1)
	REAL*4 IMAGE(IDIM,IDIM)
	COMMON/POT1/PI,QQ,PP,PASX,PASY,PASZ
C
C  Projects the density on the plane XY
C
	IF(IOPT.EQ.1) THEN
 
	DO 101 IX=0,NPTS
	DO 101 IY=0,NPTS
	RRR=0.D0
 
	DO IZ=0,NPTS
	RRR=2*RHO(IX,IY,IZ)*PASZ+RRR
	END DO
	RRR=RRR*FACT
 
	IXX=NPTS+IX+1
	IYY=NPTS+IY+1
	IMAGE(IXX,IYY)=SNGL(RRR)
	IXX=NPTS-IX+1
	IYY=IY+NPTS+1
	IMAGE(IXX,IYY)=SNGL(RRR)
	IXX=NPTS-IX+1
	IYY=NPTS-IY+1
	IMAGE(IXX,IYY)=SNGL(RRR)
	IXX=NPTS+IX+1
	IYY=NPTS-IY+1
	IMAGE(IXX,IYY)=SNGL(RRR)
 
101	CONTINUE
	ENDIF
 
C
C  Projects the density on the plane YZ
C
	IF(IOPT.EQ.2) THEN
 
	DO 102 IX=0,NPTS
	DO 102 IY=0,NPTS
	RRR=0.D0
 
	DO IZ=1,NPTS
	RRR=2*RHO(IZ,IX,IY)*PASX+RRR
	END DO
	RRR=RRR*FACT
 
	IXX=NPTS+IX+1
	IYY=NPTS+IY+1
	IMAGE(IXX,IYY)=SNGL(RRR)
	IXX=NPTS-IX+1
	IYY=IY+NPTS+1
	IMAGE(IXX,IYY)=SNGL(RRR)
	IXX=NPTS-IX+1
	IYY=NPTS-IY+1
	IMAGE(IXX,IYY)=SNGL(RRR)
	IXX=NPTS+IX+1
	IYY=NPTS-IY+1
	IMAGE(IXX,IYY)=SNGL(RRR)
 
102	CONTINUE
	ENDIF
 
C
C  Projects the density on the plane ZX
C
	IF(IOPT.EQ.3) THEN
 
	DO 103 IX=0,NPTS
	DO 103 IY=0,NPTS
	RRR=0.D0
 
	DO IZ=0,NPTS
	RRR=2*RHO(IY,IZ,IX)*PASY+RRR
	END DO
	RRR=RRR*FACT
 
	IXX=NPTS+IX+1
	IYY=NPTS+IY+1
	IMAGE(IXX,IYY)=SNGL(RRR)
	IXX=NPTS-IX+1
	IYY=IY+NPTS+1
	IMAGE(IXX,IYY)=SNGL(RRR)
	IXX=NPTS-IX+1
	IYY=NPTS-IY+1
	IMAGE(IXX,IYY)=SNGL(RRR)
	IXX=NPTS+IX+1
	IYY=NPTS-IY+1
	IMAGE(IXX,IYY)=SNGL(RRR)
 
103	CONTINUE
	ENDIF
 
	RETURN
	END
 
	SUBROUTINE PLANE(PHI,NPTS,IDIM1,IMAGE,IDIM,FACT,IOPT)
C****************************************************
C	Subroutine PLANE to create an image of the potential
C	in the plane chosen by the user.
C	 If IOPT=1 plane X,Y
C	 If IOPT=2 plane Y,Z
C	 If IOPT=3 plane Z,X
C********************************************************
	IMPLICIT DOUBLE PRECISION (A-H,O-Z)
	REAL*8 PHI(0:IDIM1,0:IDIM1,0:IDIM1)
	REAL*4 IMAGE(IDIM,IDIM)
	COMMON/POT1/PI,QQ,PP,PASX,PASY,PASZ
 
C Projection on the XY plane
C
	IF(IOPT.EQ.1)THEN
	IZ=0
	DO 101 IX=0,NPTS
	DO 101 IY=0,NPTS
	RRR=PHI(IX,IY,IZ)
	IXX=NPTS+IX+1
	IYY=NPTS+IY+1
	IMAGE(IXX,IYY)=SNGL(RRR)
	IXX=NPTS-IX+1
	IYY=IY+NPTS+1
	IMAGE(IXX,IYY)=SNGL(RRR)
	IXX=NPTS-IX+1
	IYY=NPTS-IY+1
	IMAGE(IXX,IYY)=SNGL(RRR)
	IXX=NPTS+IX+1
	IYY=NPTS-IY+1
	IMAGE(IXX,IYY)=SNGL(RRR)
101	CONTINUE
	ENDIF
 
C Projection on the YZ plane
C
	IF(IOPT.EQ.2)THEN
	IZ=0
	DO 102 IX=0,NPTS
	DO 102 IY=0,NPTS
	RRR=PHI(IZ,IX,IY)
	IXX=NPTS+IX+1
	IYY=NPTS+IY+1
	IMAGE(IXX,IYY)=SNGL(RRR)
	IXX=NPTS-IX+1
	IYY=IY+NPTS+1
	IMAGE(IXX,IYY)=SNGL(RRR)
	IXX=NPTS-IX+1
	IYY=NPTS-IY+1
	IMAGE(IXX,IYY)=SNGL(RRR)
	IXX=NPTS+IX+1
	IYY=NPTS-IY+1
	IMAGE(IXX,IYY)=SNGL(RRR)
102	CONTINUE
	ENDIF
 
C Projection on the XY plane
C
	IF(IOPT.EQ.3)THEN
	IZ=0
	DO 103 IX=0,NPTS
	DO 103 IY=0,NPTS
	RRR=PHI(IY,IZ,IX)
	IXX=NPTS+IX+1
	IYY=NPTS+IY+1
	IMAGE(IXX,IYY)=SNGL(RRR)
	IXX=NPTS-IX+1
	IYY=IY+NPTS+1
	IMAGE(IXX,IYY)=SNGL(RRR)
	IXX=NPTS-IX+1
	IYY=NPTS-IY+1
	IMAGE(IXX,IYY)=SNGL(RRR)
	IXX=NPTS+IX+1
	IYY=NPTS-IY+1
	IMAGE(IXX,IYY)=SNGL(RRR)
103	CONTINUE
	ENDIF
 
	RETURN
	END
C********************************************************************
C         IFAIL=0
C         CALL D01GAF(RAD1,Y,NPTS,S(J),ERR(J),IFAIL)
C D01GAF integrates Y(RAD1). Method of GILL and MILLER
C Comp. Journal 15, pp. 80-83, 1972. An algorithm for the integration
C of unequally spaced data.
C 	RAD1: all distinct in either ascending or descending order
C	S(J) : (output) estimate of the integral
C	ERR(J) : (output) estimate of the uncertainty in the integral
C
