Changeset 5369


Ignore:
Timestamp:
Nov 11, 2022 10:04:18 AM (5 months ago)
Author:
toby
Message:

Incorporate James Hester's updates for high angle FCJ

Location:
trunk/fsource
Files:
1 added
1 deleted
3 edited

Legend:

Unmodified
Added
Removed

  • trunk/fsource/SConstruct

    r5360 r5369  
    6868for pth in [F2PYpath,spath,os.path.normpath(os.path.join(spath,'..')),os.path.join(spath,'Scripts')]:
    6969    if not pth: continue
     70    # look for f2py3 first
    7071    if sys.platform == "win32":
    7172        program = 'f2py3.exe'
    … …  
    7677        F2PYpath,F2PYprog = os.path.split(f2pyprogram)
    7778        break
     79    # not there, try f2py
    7880    if sys.platform == "win32":
    7981        program = 'f2py.exe'
    … …  
    8486        F2PYpath,F2PYprog = os.path.split(f2pyprogram)
    8587        break
     88    # none of the above, look for f2py.py (probably obsolete)
    8689    program = 'f2py.py'
    8790    f2pyprogram = os.path.join(pth,program)
    … …  
    412415        #print ('Compile: ',file, target)
    413416        filelist.append(target)
     417    for file in glob.glob(os.path.join(sub,'*.f90')):
     418        target = env.fort(file) # connect .o files to .f90 files
     419        filelist.append(target)
    414420    #lib = Library(sub, Glob(os.path.join(sub,'*.for'))) # register library to be created
    415421    if sys.platform == "win32":

  • trunk/fsource/powsubs/psvoigt.for

    r1274 r5369  
    2626      REAL*4        GNORM               !Gaussian Normalization constant
    2727      REAL*4        COFT(6),COFN(3)
    28 
     28      REAL*4        EPS                 !Are values different
     29! Local variables saved between calls
     30      REAL*4        prev_sig,prev_gam
     31      REAL*4        eta,fwhm,frac,dsdg,dsdl,sumhm,dedf,sqsg
     32      save          eta,fwhm,frac,prev_sig,prev_gam,dsdg,dsdl,sumhm,
     33     1                   dedf,sqsg
     34     
    2935!SUBROUTINES CALLED:
    3036
    … …  
    3743      DATA COFT/1.0,2.69269,2.42843,4.47163,0.07842,1.0/
    3844      DATA COFN/1.36603,-0.47719,0.11116/
     45      DATA prev_sig/-1.0/
     46      DATA prev_gam/-1.0/
     47      DATA eps/1.0e-10/           !Threshold for recalculation
    3948
    4049!CODE:
    4150
    42       SQSG = MAX(SQRT(SIG),0.001)
    43       FWHG = STOFW*SQSG
    44       PGL = FWHG**5
    45       SUMHM = PGL
    46       DSDL = 0.0
    47       DSDG = 0.0
    48       DO ITRM=1,5
    49         PGL = PGL/FWHG
    50         DSDL = DSDL+FLOAT(ITRM)*COFT(ITRM+1)*PGL
    51         DSDG = DSDG+FLOAT(6-ITRM)*COFT(ITRM)*PGL
    52         PGL = PGL*GAM
    53         SUMHM = SUMHM+COFT(ITRM+1)*PGL
    54       END DO
    55       FWHM = EXP(0.2*LOG(SUMHM))
    56       FRAC = GAM/FWHM
    57       DEDF = 0.0
    58       PF = 1.0
    59       ETA = 0.0
    60       DO ITRM=1,3
    61         DEDF = DEDF+FLOAT(ITRM)*COFN(ITRM)*PF
    62         PF = PF*FRAC
    63         ETA = ETA+COFN(ITRM)*PF
    64       END DO
     51! Check for repeat call
     52      if (abs(prev_sig-sig) .gt. eps .or.   
     53     1  (abs(prev_gam-gam).gt.eps)) then !need to recalculate
     54         prev_sig = sig
     55         prev_gam = gam
     56         SQSG = MAX(SQRT(SIG),0.001)
     57         FWHG = STOFW*SQSG
     58         PGL = FWHG**5
     59         SUMHM = PGL
     60         DSDL = 0.0
     61         DSDG = 0.0
     62         DO ITRM=1,5
     63            PGL = PGL/FWHG
     64            DSDL = DSDL+FLOAT(ITRM)*COFT(ITRM+1)*PGL
     65            DSDG = DSDG+FLOAT(6-ITRM)*COFT(ITRM)*PGL
     66            PGL = PGL*GAM
     67            SUMHM = SUMHM+COFT(ITRM+1)*PGL
     68         END DO
     69         FWHM = EXP(0.2*LOG(SUMHM))
     70         FRAC = GAM/FWHM
     71         DEDF = 0.0
     72         PF = 1.0
     73         ETA = 0.0
     74         DO ITRM=1,3
     75            DEDF = DEDF+FLOAT(ITRM)*COFN(ITRM)*PF
     76            PF = PF*FRAC
     77            ETA = ETA+COFN(ITRM)*PF
     78         END DO
     79      end if   !end of recalculation step
    6580      CALL LORENTZ(DX,FWHM,TL,DTLDT,DTLDFW)
    6681      SIGP = (FWHM/STOFW)**2
    … …  
    8398      RETURN
    8499      END
    85 
     100     
    86101      SUBROUTINE PSVOIGT2(DX,SIG,GAM,FUNC,DFDX,DFDS,DFDG)
    87102

  • trunk/fsource/pypowder.for

    r1970 r5369  
    1       SUBROUTINE PYPSVFCJ(NPTS,DTT,TTHETA,SIG,GAM,SPH,PRFUNC)
     1      Subroutine PYPSVFCJ(NPTS,DTT,TTHETA,CDSIG,CDGAM,SPH,PRFUNC)
    22C DTT in degrees
    33C TTHETA in degrees
    4 C SPH is S/L + H/L
     4C     SPH is S/L + H/L
     5C     CDSIG Gaussian variance,centidegrees squared
     6C     CDGAM Lorenzian FWHM, centidegrees
    57C RETURNS FUNCTION ONLY
    68Cf2py intent(in) NPTS
    79Cf2py intent(in) DTT
    8 cf2py depend(NPTS) DTT
     10Cf2py depend(NPTS) DTT
    911Cf2py intent(in) TTHETA
    10 Cf2py intent(in) SIG
    11 Cf2py intent(in) GAM
     12Cf2py intent(in) CDSIG
     13Cf2py intent(in) CDGAM
    1214Cf2py intent(in) SPH
    1315Cf2py intent(out) PRFUNC
    1416Cf2py depend(NPTS) PRFUNC
    15 
    1617      REAL*4 DTT(0:NPTS-1),PRFUNC(0:NPTS-1)
    17       REAL*4 TTHETA,SIG,GAM,SPH
     18      REAL*4 TTHETA,SIG,CDSIG,CDGAM,SPH,DPRDT,DPRDG,DPRDD,DPRDLZ,DPRDS
     19      REAL*4 GAM
    1820      INTEGER*4 NPTS,I
    19       DO I=0,NPTS-1
    20         CALL PSVFCJ(DTT(I)*100.,TTHETA*100.,SIG,GAM,SPH,
    21      1    PRFUNC(I),DPRDT,SIGPART,GAMPART,SLPART)
    22       END DO
    23       RETURN
    24       END
    25 
    26       SUBROUTINE PYDPSVFCJ(NPTS,DTT,TTHETA,SIG,GAM,SPH,PRFUNC,
    27      1  DPRDT,SIGPART,GAMPART,SLPART)
     21C     CDSIG is in centidegrees squared, we must change to normal degrees
     22      SIG = CDSIG/10000.0
     23      GAM = CDGAM/100.0
     24      DO I=0,NPTS-1
     25         Call Get_Prof_Val(SIG,GAM,SPH,0.0,DTT(I)+TTHETA,TTHETA,DPRDT,
     26     1    DPRDG,DPRDD,DPRDS,DPRDLZ,PRFUNC(I))
     27          PRFUNC(I)=PRFUNC(I)/100.   !Calling code expects peak normalised in centidegrees
     28      END DO
     29      RETURN
     30      END
     31
     32      SUBROUTINE PYDPSVFCJ(NPTS,DTT,TTHETA,CDSIG,CDGAM,SPH,PRFUNC,
     33     1   DPRDT,SIGPART,GAMPART,SLPART)
    2834C DTT in degrees
    2935C TTHETA in degrees
    … …  
    3238Cf2py intent(in) NPTS
    3339Cf2py intent(in) DTT
    34 cf2py depend(NPTS) DTT
     40Cf2py depend(NPTS) DTT
    3541Cf2py intent(in) TTHETA
    36 Cf2py intent(in) SIG
    37 Cf2py intent(in) GAM
     42Cf2py intent(in) CDSIG
     43Cf2py intent(in) CDGAM
    3844Cf2py intent(in) SPH
    3945Cf2py intent(out) PRFUNC
    … …  
    4753Cf2py intent(out) SLPART
    4854Cf2py depend(NPTS) SLPART
    49 
    50       INTEGER*4 NPTS
    51       REAL*4 TTHETA,SIG,GAM,SPH
    52       REAL*4 DTT(0:NPTS-1),DPRDT(0:NPTS-1),SIGPART(0:NPTS-1),
    53      1  GAMPART(0:NPTS-1),SLPART(0:NPTS-1),PRFUNC(0:NPTS-1)
    54       DO I=0,NPTS-1
    55         CALL PSVFCJ(DTT(I)*100.,TTHETA*100.,SIG,GAM,SPH,
    56      1    PRFUNC(I),DPRDT(I),SIGPART(I),GAMPART(I),SLPART(I))
    57         DPRDT(I) = DPRDT(I)*100.
     55      INTEGER*4 NPTS,I
     56      REAL*4 TTHETA,CDSIG,SIG,CDGAM,SPH,LPART
     57      REAL*4 GAM
     58      REAL*4 DTT(0:NPTS-1),DPRDT(0:NPTS-1),SIGPART(0:NPTS-1)
     59      REAL*4 GAMPART(0:NPTS-1),SLPART(0:NPTS-1),PRFUNC(0:NPTS-1)
     60      SIG = CDSIG/10000.
     61      GAM = CDGAM/100.
     62      DO I=0,NPTS-1
     63         Call Get_Prof_Val(SIG,GAM,SPH,0.0,DTT(I)+TTHETA,TTHETA,
     64     1      DPRDT(I),SIGPART(I),GAMPART(I),SLPART(I),LPART, PRFUNC(I))
     65         ! Calling code expects all values to be for a peak normalised in centidegrees
     66         SIGPART(I)=SIGPART(I)/1.0e6
     67         GAMPART(I)=GAMPART(I)/1.0e4
     68         SLPART(I)=SLPART(I)/100.
     69         PRFUNC(I)=PRFUNC(I)/100.
     70         DPRDT(I)=DPRDT(I)/100.
    5871      END DO
    5972      RETURN
    … …  
    6578Cf2py intent(in) NPTS
    6679Cf2py intent(in) DTT
    67 cf2py depend(NPTS) DTT
     80Cf2py depend(NPTS) DTT
    6881Cf2py intent(in) SIG
    6982Cf2py intent(in) GAM
    … …  
    7689      DO I=0,NPTS-1
    7790        CALL PSVOIGT(DTT(I)*100.,SIG,GAM,
    78      1    PRFUNC(I),DPRDT,SIGPART,GAMPART)
     91     1     PRFUNC(I),DPRDT,SIGPART,GAMPART)
    7992      END DO
    8093      RETURN
    … …  
    8295
    8396      SUBROUTINE PYDPSVOIGT(NPTS,DTT,SIG,GAM,PRFUNC,
    84      1  DPRDT,SIGPART,GAMPART)
     97     1   DPRDT,SIGPART,GAMPART)
    8598C DTT in degrees
    8699C RETURNS FUNCTION & DERIVATIVES
    87100Cf2py intent(in) NPTS
    88101Cf2py intent(in) DTT
    89 cf2py depend(NPTS) DTT
     102Cf2py depend(NPTS) DTT
    90103Cf2py intent(in) SIG
    91104Cf2py intent(in) GAM
    … …  
    102115      REAL*4 SIG,GAM
    103116      REAL*4 DTT(0:NPTS-1),DPRDT(0:NPTS-1),SIGPART(0:NPTS-1),
    104      1  GAMPART(0:NPTS-1),PRFUNC(0:NPTS-1)
     117     1   GAMPART(0:NPTS-1),PRFUNC(0:NPTS-1)
    105118      DO I=0,NPTS-1
    106119        CALL PSVOIGT(DTT(I)*100.,SIG,GAM,
    107      1    PRFUNC(I),DPRDT(I),SIGPART(I),GAMPART(I))
     120     1     PRFUNC(I),DPRDT(I),SIGPART(I),GAMPART(I))
    108121        DPRDT(I) = DPRDT(I)*100.
    109122      END DO
    … …  
    118131Cf2py intent(in) NPTS
    119132Cf2py intent(in) DTT
    120 cf2py depend(NPTS) DTT
     133Cf2py depend(NPTS) DTT
    121134Cf2py intent(in) TTHETA
    122135Cf2py intent(in) SIG
    … …  
    131144      DO I=0,NPTS-1
    132145        CALL PSVFCJO(DTT(I)*100.,TTHETA*100.,SIG,GAM,SPH/2.0,SPH/2.0,
    133      1    PRFUNC(I),DPRDT,SIGPART,GAMPART,SLPART,HLPART)
     146     1     PRFUNC(I),DPRDT,SIGPART,GAMPART,SLPART,HLPART)
    134147      END DO
    135148      RETURN
    … …  
    144157Cf2py intent(in) NPTS
    145158Cf2py intent(in) DTT
    146 cf2py depend(NPTS) DTT
     159Cf2py depend(NPTS) DTT
    147160Cf2py intent(in) TTHETA
    148161Cf2py intent(in) SIG
    … …  
    163176      REAL*4 TTHETA,SIG,GAM,SHL
    164177      REAL*4 DTT(0:NPTS-1),DPRDT(0:NPTS-1),SIGPART(0:NPTS-1),
    165      1  GAMPART(0:NPTS-1),SLPART(0:NPTS-1),PRFUNC(0:NPTS-1)
     178     1   GAMPART(0:NPTS-1),SLPART(0:NPTS-1),PRFUNC(0:NPTS-1)
    166179      DO I=0,NPTS-1
    167180        CALL PSVFCJO(DTT(I)*100.,TTHETA*100.,SIG,GAM,SHL/2.,SHL/2.,
    168      1    PRFUNC(I),DPRDT(I),SIGPART(I),GAMPART(I),SPART,HPART)
     181     1     PRFUNC(I),DPRDT(I),SIGPART(I),GAMPART(I),SPART,HPART)
    169182          SLPART(I) = SPART
    170183        DPRDT(I) = DPRDT(I)*100.
    … …  
    178191Cf2py intent(in) NPTS
    179192Cf2py intent(in) DTT
    180 cf2py depend(NPTS) DTT
     193Cf2py depend(NPTS) DTT
    181194Cf2py intent(in) ALP
    182195Cf2py intent(in) BET
    … …  
    189202      REAL*4 ALP,BET,SIG,GAM,SHL
    190203      REAL*4 DTT(0:NPTS-1),PRFUNC(0:NPTS-1),DPRDT,ALPPART,
    191      1  BETPART,SIGPART,GAMPART
     204     1   BETPART,SIGPART,GAMPART
    192205      DO I=0,NPTS-1
    193206        CALL EPSVOIGT(DTT(I),ALP,BET,SIG,GAM,PRFUNC(I),DPRDT,
    194      1    ALPPART,BETPART,SIGPART,GAMPART)
     207     1     ALPPART,BETPART,SIGPART,GAMPART)
    195208      END DO
    196209      RETURN
    … …  
    198211
    199212      SUBROUTINE PYDEPSVOIGT(NPTS,DTT,ALP,BET,SIG,GAM,PRFUNC,
    200      1  DPRDT,ALPPART,BETPART,SIGPART,GAMPART)
     213     1   DPRDT,ALPPART,BETPART,SIGPART,GAMPART)
    201214C DTT in microsec
    202215C RETURNS FUNCTION & DERIVATIVES
    203216Cf2py intent(in) NPTS
    204217Cf2py intent(in) DTT
    205 cf2py depend(NPTS) DTT
     218Cf2py depend(NPTS) DTT
    206219Cf2py intent(in) ALP
    207220Cf2py intent(in) BET
    … …  
    224237      REAL*4 ALP,BET,SIG,GAM,SHL
    225238      REAL*4 DTT(0:NPTS-1),DPRDT(0:NPTS-1),ALPPART(0:NPTS-1),
    226      1  BETPART(0:NPTS-1),SIGPART(0:NPTS-1),
    227      1  GAMPART(0:NPTS-1),PRFUNC(0:NPTS-1)
     239     1   BETPART(0:NPTS-1),SIGPART(0:NPTS-1),
     240     1   GAMPART(0:NPTS-1),PRFUNC(0:NPTS-1)
    228241      DO I=0,NPTS-1
    229242        CALL EPSVOIGT(DTT(I),ALP,BET,SIG,GAM,PRFUNC(I),DPRDT(I),
    230      1    ALPPART(I),BETPART(I),SIGPART(I),GAMPART(I))
     243     1     ALPPART(I),BETPART(I),SIGPART(I),GAMPART(I))
    231244      END DO
    232245      RETURN
    … …  
    261274Cf2py intent(in) NIN
    262275Cf2py intent(in)  XIN
    263 cf2py depend(NIN) XIN
     276Cf2py depend(NIN) XIN
    264277Cf2py intent(in)  YIN
    265 cf2py depend(NIN) YIN
     278Cf2py depend(NIN) YIN
    266279Cf2py intent(in) NOUT
    267280Cf2py intent(in)   XOUT
    268 cf2py depend(NOUT) XOUT
     281Cf2py depend(NOUT) XOUT
    269282Cf2py intent(out)  YOUT
    270 cf2py depend(NOUT) YOUT
     283Cf2py depend(NOUT) YOUT
    271284
    272285      INTEGER NIN,NOUT
Note: See TracChangeset for help on using the changeset viewer.