      PROGRAM SRPLOT                                                    
*----------------------------------------------------------------------*
*                                                                      *
*    Written by A. Fasso', Sept. 30, 1993                              *
*    SRPLOT calculates a differential or an integral synchrotron       *
*    radiation spectrum and prepares a TOPDRAW data file to display    *
*    it. Options are photon or power spectra. The user may input the   *
*    critical energy, if known. Otherwise, the program calculates it   *
*    from the bending radius and from the electron beam energy.        *
*    The output file must still be edited to change the Y limits, but  *
*    the program provides the max. and min. Y to guide the choice.     *
*    The algorithm used is that of Umstaetter (CERN), as implemented   *
*    by G.R. Stevenson                                                 *
*                                                                      *
*----------------------------------------------------------------------*
      IMPLICIT DOUBLE PRECISION (A-H, O-Z)                              
      EXTERNAL SYNRAD                                                   
      CHARACTER*80 TERM, FILOUT*20                                      
      LOGICAL INTGRL, NORMAL                                            
      PARAMETER(NUMPOI = 100) ! 100 pts per decade for smooth spect.    
      PARAMETER(NUMDEC = 9)    ! 9 decades should suffice               
      DOUBLE PRECISION EARRAY(NUMPOI*NUMDEC), YARRAY(NUMPOI*NUMDEC)     
      DATA INTGRL/.FALSE./                                              
*--------------------  Questions and answers  -------------------------*
      WRITE(*,*) '           ***** Answer Q or q to quit *****'         
      WRITE(*,*) '           Give the name of the output file '         
      READ(*,'(A80)') TERM                                              
      IF(INDEX(TERM,'Q') .GT. 0 .OR. INDEX(TERM,'q') .GT. 0) STOP       
      READ(TERM,'(A20)') FILOUT                                         
      OPEN ( UNIT=1, STATUS='UNKNOWN', FILE=FILOUT)                
*-----------------                                                      
    1 CONTINUE                                                          
      WRITE(*,*) '         Do you want a photon energy spectrum (1),'   
      WRITE(*,*) '          a photon lethargy (logE) spectrum (2),'     
      WRITE(*,*) '               a power energy spectrum (3),'          
      WRITE(*,*) '            or a power lethargy spectrum (4)?'        
      WRITE(*,*) '     (if you want an integral spectrum select 1 or 3)'
      READ(*,'(A80)') TERM                                              
      IF(INDEX(TERM,'Q') .GT. 0 .OR. INDEX(TERM,'q') .GT. 0) STOP       
      READ(TERM,*) IANSW1                                               
*-----------------                                                      
      IF(IANSW1 .EQ. 1 .OR. IANSW1 .EQ. 3) THEN                         
        WRITE(*,*)                                                      
     &   'Do you want a differential (1) or an integral spectrum (2) ?' 
        READ(*,'(A80)') TERM                                            
        IF(INDEX(TERM,'Q') .GT. 0 .OR. INDEX(TERM,'q') .GT. 0) STOP     
        READ(TERM,*) IANSW3                                             
*       -----------                                                     
        IF(IANSW3 .EQ. 2) THEN                                          
          INTGRL =.TRUE.                                                
          WRITE(*,*) '         Do you want it normalized? (Y/N)'        
          READ(*,'(A80)') TERM                                          
          IF(INDEX(TERM,'Q') .GT. 0 .OR. INDEX(TERM,'q') .GT. 0) STOP   
          IF(INDEX(TERM,'Y') .GT. 0 .OR. INDEX(TERM,'y') .GT. 0) THEN       
            NORMAL = .TRUE.                                             
          ELSE
            NORMAL = .FALSE.
          ENDIF
        ENDIF                                                           
      ENDIF                                                             
*-----------------                                                      
      WRITE(*,*)'      Give the lower energy limit of the plot (in MeV)'
      READ(*,'(A80)') TERM                                              
      IF(INDEX(TERM,'Q') .GT. 0 .OR. INDEX(TERM,'q') .GT. 0) STOP       
      READ(TERM,*) EMIN                                                 
      WRITE(*,*)'     Give the higher energy limit of the plot (in MeV)'
      READ(*,'(A80)') TERM                                              
      IF(INDEX(TERM,'Q') .GT. 0 .OR. INDEX(TERM,'q') .GT. 0) STOP       
      READ(TERM,*) EMAX                                                 
*-----------------              
      WRITE(1,*) 'set font duplex'                                  
      WRITE(1,*) 'set title size 1.5'
      IF(IANSW1 .EQ. 1) THEN                                            
        IF(INTGRL) THEN                                                 
          WRITE(1,*) 'title top ',                                      
     &    '''Synchrotron radiation integral photon spectrum'''          
        ELSE                                                            
        WRITE(1,*) 'title top ',                                        
     &    '''Synchrotron radiation photon spectrum'''                   
        ENDIF                                                           
        WRITE(1,*) 'title bottom ''Energy (MeV)'''                      
        IF(INTGRL) THEN                                                 
          WRITE(1,*) 'title left ',                                     
     &    '''                  photon fraction below given energy '''   
        ELSE                                                            
          WRITE(1,*) 'title left ',                                     
     &    '''                  dN/dk (photons per m per MeV)'''         
        ENDIF                                                           
        WRITE(1,*) 'set scale x log y log'                              
        WRITE(1,'(1X,''set limits x '',1P,G12.4,'' to '',G12.4,         
     &  '' y 1.E-30 to 1.5'' )') EMIN, EMAX                             
      ELSEIF(IANSW1 .EQ. 2) THEN                                        
        WRITE(1,*) 'title top ''Synchrotron radiation photon spectrum'''
        WRITE(1,*) 'title bottom ''Energy (MeV)'''                      
        WRITE(1,*) 'title left ',                                       
     &  '''              E * dN/dk (photons per m per unit lethargy)''' 
        WRITE(1,*) 'set scale x log y lin'                              
        WRITE(1,'(1X,''set limits x '',1P,G12.4,'' to '',G12.4,         
     &  '' y 0.0 to 1.5'' )') EMIN, EMAX                                
      ELSEIF(IANSW1 .EQ. 3) THEN                                        
        IF(INTGRL) THEN                                                 
          WRITE(1,*) 'title top ',                                      
     &    '''Synchrotron radiation integral power spectrum'''           
        ELSE                                                            
          WRITE(1,*) 'title top ',                                      
     &    '''Synchrotron radiation power spectrum'''                    
        ENDIF                                                           
        WRITE(1,*) 'title bottom ''Energy (MeV)'''                      
        IF(INTGRL) THEN                                                 
          WRITE(1,*) 'title left ',                                     
     &    '''                  power fraction below given energy'''     
        ELSE                                                            
          WRITE(1,*) 'title left ',                                     
     &    '''                  dE/dk (MeV per m per MeV)'''             
        ENDIF                                                           
        WRITE(1,*) 'set scale x log y log'                              
        WRITE(1,'(1X,''set limits x '',1P,G12.4,'' to '',G12.4,         
     &  '' y 1.E-30 to 1.5'' )') EMIN, EMAX                             
      ELSEIF(IANSW1 .EQ. 4) THEN                                        
        WRITE(1,*) 'title top ''Synchrotron radiation power spectrum''' 
        WRITE(1,*) 'title bottom ''Energy (MeV)'''                      
        WRITE(1,*) 'title left ',                                       
     &  '''              E * dE/dk (MeV per m per unit lethargy)'''     
        WRITE(1,*) 'set scale x log y lin'                              
        WRITE(1,'(1X,''set limits x '',1P,G12.4,'' to '',G12.4,         
     &  '' y 0.0 to 1.5'' )') EMIN, EMAX                                
      ELSE                                                              
        WRITE(*,*)                                                      
     &   'The only allowed answers are 1, 2, 3, 4 or Q. Try again'      
        GO TO 1                                                         
      ENDIF                                                             
      WRITE(1,*) 'set order x y'                                        
*-----------------                                                      
      WRITE(*,*) '           What is the electron beam energy in GeV?'  
      READ(*,'(A80)') TERM                                              
      IF(INDEX(TERM,'Q') .GT. 0 .OR. INDEX(TERM,'q') .GT. 0) STOP       
      READ(TERM,*) BEAMEN                                               
    2 CONTINUE                                                          
      WRITE(*,*) '          Do you already know the critical energy (1)'
      WRITE(*,*) '            or do you want me to calculate it (2) ?'  
      READ(*,'(A80)') TERM                                              
      IF(INDEX(TERM,'Q') .GT. 0 .OR. INDEX(TERM,'q') .GT. 0) STOP       
      READ(TERM,*) IANSW2                                               
      IF(IANSW2 .EQ. 1) THEN                                            
        WRITE(*,*) '       Give the value of the critical energy in keV'
        READ(*,'(A80)') TERM                                            
        IF(INDEX(TERM,'Q') .GT. 0 .OR. INDEX(TERM,'q') .GT. 0) STOP     
        READ(TERM,*) ECRIT                                              
        ECRIT = ECRIT * 1.D-3     !   Internally we use MeV             
      ELSEIF(IANSW2 .EQ. 2) THEN                                        
        WRITE(*,*) '           Give the bending radius in meters'       
        READ(*,'(A80)') TERM                                            
        IF(INDEX(TERM,'Q') .GT. 0 .OR. INDEX(TERM,'q') .GT. 0) STOP     
        READ(TERM,*) RADIUS                                             
        ECRIT=2.218D-3*BEAMEN**3/RADIUS                                 
      ELSE                                                              
        WRITE(*,*) '  The only allowed answers are 1, 2 or Q. Try again'
        GO TO 2                                                         
      ENDIF                                                             
*----------------------------------------------------------------------*
*   Calculates the number of synchrotron radiation photons per metre   *
*   of orbit per electron per energy bin (MeV units)                   *
*----------------------------------------------------------------------*
      FACTOR = 1.775D3/(BEAMEN*BEAMEN)                                  
      ENERGY = EMIN                                                     
      EN1 = ENERGY/ECRIT                                                
      IF(IANSW1 .EQ. 1) THEN                                            
        Y1=SYNRAD(EN1)*FACTOR                  !  Photon energy spectrum
      ELSEIF(IANSW1 .EQ. 2 .OR. IANSW1 .EQ. 3) THEN                     
*                               Photon lethargy or power energy spectrum
        Y1=SYNRAD(EN1)*FACTOR*ENERGY                                    
      ELSEIF(IANSW1 .EQ. 4) THEN                                        
        Y1=SYNRAD(EN1)*FACTOR*ENERGY**2       !  Power lethargy spectrum
      ENDIF                                                             
      IF((IANSW1 .EQ. 2 .OR. IANSW1 .EQ. 4                              
     &   .OR. Y1 .GT. 1.D-40) .AND. .NOT. INTGRL ) THEN                 
        WRITE(1,300) ENERGY, Y1                                         
 300    FORMAT(1P,4E15.5)                                               
        YMIN = Y1                                                       
      ELSE                                                              
        YMIN = 1.D-40                                                   
      ENDIF                                                             
      YMAX = Y1                                                         
      IF(INTGRL) YSUM = 0.D0                                            
      IDCMIN = NINT(LOG10(EMIN))                                        
      IDCMAX = NINT(LOG10(EMAX))                                        
      RATIO = 10.D0**(1.D0/DBLE(NUMPOI))                                
      KOUNT = 0                                                         
      DO IDECAD = IDCMIN, IDCMAX                                        
*       DECFAC = 10.D0**IDECAD                                          
        DO J = 1, NUMPOI                                                
*         ENEW = ENERGY + DECFAC / DBLE(NUMPOI)                         
          ENEW = ENERGY * RATIO                                         
          IF(ENEW .GT. EMAX) GO TO 3                                    
          KOUNT = KOUNT + 1                                             
          IF(INTGRL) DELTA = ENEW - ENERGY                              
          ENERGY = ENEW                                                 
          EN1 = ENERGY/ECRIT                                            
          IF(IANSW1 .EQ. 1) THEN                                        
            Y1=SYNRAD(EN1)*FACTOR              !  Photon energy spectrum
          ELSEIF(IANSW1 .EQ. 2 .OR. IANSW1 .EQ. 3) THEN                 
*                               Photon lethargy or power energy spectrum
            Y1=SYNRAD(EN1)*FACTOR*ENERGY                                
          ELSEIF(IANSW1 .EQ. 4) THEN                                    
            Y1=SYNRAD(EN1)*FACTOR*ENERGY**2   !  Power lethargy spectrum
          ENDIF                                                         
          IF(INTGRL) THEN                                               
            YSUM = YSUM + DELTA * Y1                                    
            EARRAY(KOUNT) = ENERGY                                      
            YARRAY(KOUNT) = YSUM                                        
            IF(IDECAD .EQ. IDCMIN .AND. J .EQ. 1) YSUMIN = YSUM         
          ELSEIF(IANSW1 .EQ. 2 .OR. IANSW1 .EQ. 4                       
     &                                         .OR. Y1 .GT. 1.D-40) THEN
            WRITE(1,300) ENERGY, Y1                                     
            IF(Y1 .LT. YMIN) YMIN = Y1                                  
          ENDIF                                                         
          IF(Y1 .GT. YMAX) YMAX = Y1                                    
        END DO                                                          
      END DO                                                            
    3 CONTINUE                                                          
      IF(INTGRL) THEN                                                   
        DO K = 1, KOUNT                                                 
          IF(NORMAL) YARRAY(K) = YARRAY(K) / YSUM                       
          IF(YARRAY(K) .GT. 1.D-40)                                     
     &      WRITE(1,'(1X,1P,G28.20,G28.20)') EARRAY(K), YARRAY(K)       
        END DO                                                          
      ENDIF                                                             
      WRITE(1,*) 'join'                                                 
      IF(INTGRL) THEN                                                   
        YLO = YSUMIN                                                    
        YHI = YSUM                                                      
      ELSE                                                              
        YLO = YMIN                                                      
        YHI = YMAX                                                      
      ENDIF                                                             
*     WRITE(1,400) EMIN*10.D0, YHI*0.75D0, ECRIT*1.D3                   
* 400 FORMAT('TITLE ',1P,2G15.5,'DATA ''E0c1 =',G11.4,' keV'' ')        
*     WRITE(1,500)                                                      
* 500 FORMAT(                   'CASE ''UXLX                 ''')       
      CLOSE(UNIT = 1)                                                   
      WRITE(*,*) 'The TOPDRAW file is ready, ',                         
     &     'but you may need to change the Y limits:'                   
      IF(IANSW1 .EQ. 1 .OR. IANSW1 .EQ. 3)                              
     &  WRITE(*,'(''           The minimum Y-value is '',1P,G12.4)') YLO
      IF(.NOT. NORMAL)                                                  
     &  WRITE(*,'(''           The maximum Y-value is '',1P,G12.4)') YHI
      END                                                               
      DOUBLE PRECISION FUNCTION SYNRAD(X)                               
      DOUBLE PRECISION A,B,P,Q,X,Y,Z                                    
      RX=X                                                              
      IF(RX.LT.6.) GO TO 1                                              
      IF(RX.GT.174.) GO TO 2                                            
      Z=20.D0/X-2.D0                                                    
      A=               +.00000000000000000001D0                         
      B=     Z*A       -.00000000000000000002D0                         
      A=     Z*B-A     +.00000000000000000006D0                         
      B=     Z*A-B     -.00000000000000000020D0                         
      A=     Z*B-A     +.00000000000000000066D0                         
      B=     Z*A-B     -.00000000000000000216D0                         
      A=     Z*B-A     +.00000000000000000721D0                         
      B=     Z*A-B     -.00000000000000002443D0                         
      A=     Z*B-A     +.00000000000000008441D0                         
      B=     Z*A-B     -.00000000000000029752D0                         
      A=     Z*B-A     +.00000000000000107116D0                         
      B=     Z*A-B     -.00000000000000394564D0                         
      A=     Z*B-A     +.00000000000001489474D0                         
      B=     Z*A-B     -.00000000000005773537D0                         
      A=     Z*B-A     +.00000000000023030657D0                         
      B=     Z*A-B     -.00000000000094784973D0                         
      A=     Z*B-A     +.00000000000403683207D0                         
      B=     Z*A-B     -.00000000001785432348D0                         
      A=     Z*B-A     +.00000000008235329314D0                         
      B=     Z*A-B     -.00000000039817923621D0                         
      A=     Z*B-A     +.00000000203088939238D0                         
      B=     Z*A-B     -.00000001101482369622D0                         
      A=     Z*B-A     +.00000006418902302372D0                         
      B=     Z*A-B     -.00000040756144386809D0                         
      A=     Z*B-A     +.00000287536465397527D0                         
      B=     Z*A-B     -.00002321251614543524D0                         
      A=     Z*B-A     +.00022505317277986004D0                         
      B=     Z*A-B     -.00287636803664026799D0                         
      A=     Z*B-A     +.06239591359332750793D0                         
      P=.5D0*Z*A-B    +1.06552390798340693166D0                         
      SYNRAD=P*DSQRT(1.57079 63267 94896 61192 3132169164D0/X)/DEXP(X)  
      RETURN                                                            
    1 Z=X**2/16.D0-2.D0                                                 
      B=                .00000000000000000012D0                         
      A=     Z*B  +     .00000000000000000460D0                         
      B=     Z*A-B+     .00000000000000031738D0                         
      A=     Z*B-A+     .00000000000002004426D0                         
      B=     Z*A-B+     .00000000000111455474D0                         
      A=     Z*B-A+     .00000000005407460944D0                         
      B=     Z*A-B+     .00000000226722011790D0                         
      A=     Z*B-A+     .00000008125130371644D0                         
      B=     Z*A-B+     .00000245751373955212D0                         
      A=     Z*B-A+     .00006181256113829740D0                         
      B=     Z*A-B+     .00127066381953661690D0                         
      A=     Z*B-A+     .02091216799114667278D0                         
      B=     Z*A-B+     .26880346058164526514D0                         
      A=     Z*B-A+    2.61902183794862213818D0                         
      B=     Z*A-B+   18.65250896865416256398D0                         
      A=     Z*B-A+   92.95232665922707542088D0                         
      B=     Z*A-B+  308.15919413131586030542D0                         
      A=     Z*B-A+  644.86979658236221700714D0                         
      P=.5D0*Z*A-B+  414.56543648832546975110D0                         
      A=                .00000000000000000004D0                         
      B=     Z*A  +     .00000000000000000289D0                         
      A=     Z*B-A+     .00000000000000019786D0                         
      B=     Z*A-B+     .00000000000001196168D0                         
      A=     Z*B-A+     .00000000000063427729D0                         
      B=     Z*A-B+     .00000000002923635681D0                         
      A=     Z*B-A+     .00000000115951672806D0                         
      B=     Z*A-B+     .00000003910314748244D0                         
      A=     Z*B-A+     .00000110599584794379D0                         
      B=     Z*A-B+     .00002581451439721298D0                         
      A=     Z*B-A+     .00048768692916240683D0                         
      B=     Z*A-B+     .00728456195503504923D0                         
      A=     Z*B-A+     .08357935463720537773D0                         
      B=     Z*A-B+     .71031361199218887514D0                         
      A=     Z*B-A+    4.26780261265492264837D0                         
      B=     Z*A-B+   17.05540785795221885751D0                         
      A=     Z*B-A+   41.83903486779678800040D0                         
      Q=.5D0*Z*A-B+   28.41787374362784178164D0                         
      Y=X**6.66666 66666 66666 66666 66666 66666D-1                     
      SYNRAD=(P/Y-Q*Y-1.D0)*1.81379 93642 34217 84215 53078 8143D0      
      RETURN                                                            
    2 SYNRAD=0.D0                                                       
      RETURN                                                            
      END                                                               
