C @(#)franz2d.for 17.1.1.1 (ES0-DMD) 01/25/02 17:10:43 C=========================================================================== C Copyright (C) 1995 European Southern Observatory (ESO) C C This program is free software; you can redistribute it and/or C modify it under the terms of the GNU General Public License as C published by the Free Software Foundation; either version 2 of C the License, or (at your option) any later version. C C This program is distributed in the hope that it will be useful, C but WITHOUT ANY WARRANTY; without even the implied warranty of C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the C GNU General Public License for more details. C C You should have received a copy of the GNU General Public C License along with this program; if not, write to the Free C Software Foundation, Inc., 675 Massachusetss Ave, Cambridge, C MA 02139, USA. C C Corresponding concerning ESO-MIDAS should be addressed as follows: C Internet e-mail: midas@eso.org C Postal address: European Southern Observatory C Data Management Division C Karl-Schwarzschild-Strasse 2 C D 85748 Garching bei Muenchen C GERMANY C=========================================================================== C SUBROUTINE FRANZ2(X,NP,PARAM,Y1,DERIV) C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C.COPYRIGHT: Copyright (c) 1987 European Southern Observatory, C all rights reserved C C.VERSION: 1.0 ESO-FORTRAN Conversion, AA 18:01 - 21 DEC 1987 C C.LANGUAGE: F77+ESOext C C.AUTHOR: J.D.PONZ C C-------------------------------------------------------------- C********************************************** C C IMPLICIT NONE C 2-DIM ROTATIONALLY SYMMETRIC FRANZ-FUNCTION C C********************************************** C C Author: O.-G. Richter, ESO Garching C INTEGER NP REAL X(2) DOUBLE PRECISION Y1,PARAM(NP),DERIV(NP) C DOUBLE PRECISION A,A1,A2,B,C,D,E,F,G,H,O,P,Q C ! 1 / Sigma1 C = 1.0D0/ (DABS(PARAM(4))+1.0D-14) ! 1 / Sigma2 D = 1.0D0/PARAM(6) ! x - xc A1 = DBLE(X(1)) - PARAM(2) ! (y - yc) / alfa A2 = (DBLE(X(2))-PARAM(3))/PARAM(7) O = A1*A2 B = A1*A1 + A2*A2 - 2.0D0*PARAM(8)*O IF (B.LT.0.0D0) B = 0.0D0 ! Distance from center B = DSQRT(B) E = 1.0D0 + B*D ! Exponent F = E*PARAM(5) ! Avoid singular value G = B*C + 1.0D-34 ! \ H = -DLOG(G)*F ! | IF (H.GE.+80.0D0) H = 1.0D37 ! | IF (DABS(H).LT.80.0D0) H = DEXP(H) ! | IF (H.LE.-80.0D0) H = 0.0D0 ! | F = 0.0D0 ! | Computation of Q = 0.0D0 ! > the denominator IF (B.EQ.0.0D0) GO TO 10 ! | avoiding overflow F = PARAM(6)/B ! | Q = DLOG(G) ! | 10 P = 0.0D0 ! | IF (H.EQ.0.0D0) GO TO 20 ! / H = 1.0D0/H ! Denominator P = 1.0D0/ (1.0D0+H) 20 DERIV(1) = P ! Function value Y1 = PARAM(1)*P P = Y1*P*H DERIV(5) = -P*E*Q P = P*PARAM(5) DERIV(4) = P*E*C P = P*D A = -P* (PARAM(4)* (1.0D0+F)+Q)/ (B+1.0D-35) DERIV(2) = A* (PARAM(8)*A2-A1) DERIV(3) = A* (PARAM(8)*A1-A2)/PARAM(7) DERIV(6) = P*B*Q*D DERIV(7) = A*A2* (PARAM(8)*A1-A2)/PARAM(7) DERIV(8) = -A*O C RETURN END