PROGRAM DEPROJELL
C
C   Program to obtain analytically deprojected estimates for the semi-major
C   axis, ellipticity, and position angle of an ellipse seen in projection,
C   following Gadotti et al. (2007, MNRAS, in press -- see their Appendix A).
C   The user should provide the name of a file (up to 25 letters or symbols)
C   which contains, for each line, four numbers separated by spaces, in the
C   following order: projected, measured semi-major axis and ellipticity
C   of the ellipse (usually the one defining a galactic bar), the difference
C   between the position angle of the ellipse and that of the line of nodes
C   (in degrees) and the inclination angle also in degrees (zero meaning no
C   inclination).
C
C   The code returns for each projected ellipse at each line of the input
C   file the deprojected values for the semi-major axis, ellipticity and
C   the angle between the ellipse and the line of nodes. The semi-major
C   axis output is in the same units as the input. The position angle output
C   is in degrees (as the input).
C
C   Although this code has been thoroughly tested you should use it
C   carefully (as any other code, by the way!). I give no warranties
C   whatsoever! If you are using it, or find mistakes, please notify me at
C
C                 ---> dimitri@mpa-garching.mpg.de <---
C
C   If you want to improve it, go ahead, but document your changes and
C   let me know!
C
C   Finally, please, if you publish any results using this code you
C
C                   *** MUST CITE ITS MAIN PAPER ***
C
C          *** Gadotti, D. A., et al. 2007, MNRAS, 381, 943 ***
C
C
C
C                                                  Dimitri Gadotti
C                                                  Munich, June, 2007
C
      IMPLICIT NONE
      CHARACTER*25 fname ! fname is the name of the user input file
      REAL smaj0 ! projected, measured semi-major axis of ellipse
      REAL ell0 ! projected, measured ellipticity
      REAL smaj ! deprojected semi-major axis of ellipse
      REAL smin ! deprojected semi-minor axis
      REAL ell ! deprojected ellipticity
      REAL pa ! deprojected position angle
      REAL alpha ! alpha
      REAL i ! i
      REAL Ap ! A'
      REAL Bp ! B'
      REAL Cp ! C'
      REAL s1 ! s1
      REAL s2 ! s2
      REAL pi ! Pi
      REAL a,b ! to intermediate estimation of smaj and ell
      REAL sint ! to intermediate estimation of s1 and s2

      WRITE (*,*) 'Please give input file name:'
      READ (*,'(1A25)') fname
      OPEN(UNIT=1,FILE=fname,STATUS='OLD')
10    READ(1,*,END=20) smaj0,ell0,alpha,i
      pi=4.0*ATAN(1.0) ! define pi
      alpha=alpha*pi/180. ! transform degrees to radians
      i=i*pi/180. ! transform degrees to radians
      a=smaj0
      b=(1.-ell0)*a
      Ap=COS(alpha)**2./a**2.+SIN(alpha)**2./b**2.
      Bp=COS(alpha)*SIN(alpha)/a**2.-COS(alpha)*SIN(alpha)/b**2.
      Bp=Bp*COS(i)
      Cp=SIN(alpha)**2./a**2.+COS(alpha)**2./b**2.
      Cp=Cp*COS(i)**2.
      sint=SQRT(1.+4.*Bp**2./(Ap-Cp)**2.)
      s1=2.*(Ap*Cp-Bp**2.)/((Bp**2.-Ap*Cp)*((Cp-Ap)*sint-Cp-Ap))
      s1=SQRT(s1)
      s2=2.*(Ap*Cp-Bp**2.)/((Bp**2.-Ap*Cp)*((Ap-Cp)*sint-Cp-Ap))
      s2=SQRT(s2)
      smaj=MAX(s1,s2)
      smin=MIN(s1,s2)
      ell=1.-smin/smaj
      pa=-0.5*ATAN(2.*Bp/(Cp-Ap))
      IF(smaj.EQ.s2) pa=pi/2.+pa
      pa=pa*180./pi
      WRITE(*,*) smaj,ell,pa
      GOTO 10
20    WRITE (*,*) 'Finished successfully!'
      CLOSE (1)
      END