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