C                                                                      
C DOPPLER SHIFT                       
C                                                                       
C NON-RELATIVISTIC KINEMATIC                                            
C 
      common/input/A1,A2,EP,THGA,PHIGA,VCM,THETA,PHI
      common/output/THCM,V1,V2,COSPG1,COSPG2,COSTG1,COSTG2,EGAMMA
C
      PI=3.1415927                                                      
C
C
C satangd output >>> OUTPUT.GRF <<<
C check the output file!!!
C
      open (21,file='OUTPUT1.GRF')
      open (22,file='OUTPUT2.GRF')
      open (23,file='OUTPUT3.GRF')
      open (24,file='OUTPUT4.GRF')
      write(21,200)
      write(22,200)
      write(23,200)
      write(24,200)
  200 format(3H' '/20HX:th^l$a$b$ (degree)/
C     128HY:Eg$a$m$m$a$/Eg$a$m$m$a$0$ )
     119HY:Eg$a$m$m$a$(MeV) )
      write(21,201)
  201 format(2Hc:/42Hc: projectile measured, projectile excited
     1/2Hc:)
      write(22,202)
  202 format(2Hc:/46Hc: projectile measured, target nucleus excited
     1/2Hc:)
      write(23,203)
  203 format(2Hc:/50Hc: target nucleus measured, target nucleus excited
     1/2Hc:)
      write(24,204)
  204 format(2Hc:/46Hc: target nucleus measured, projectile excited
     1/2Hc:)
      write(21,205)
  205 format(39HH:    X      Y,bl0 Y1,bd0 Y2,bm0 Y3,bp0)
      write(22,206)
  206 format(39HH:    X      Y,bl0 Y1,bd0 Y2,bm0 Y3,bp0)
      write(23,207)
  207 format(39HH:    X      Y,rl0 Y1,rd0 Y2,rm0 Y3,rp0)
      write(24,208)
  208 format(39HH:    X      Y,rl0 Y1,rd0 Y2,rm0 Y3,rp0)
C
      WRITE(6,100)                                                      
 100  FORMAT(1X,' A1 = ')                                              
      READ(5,*) A1                                                      
      WRITE(6,101)                                                      
 101  FORMAT(1X,' A2 = ')                                               
      READ(5,*) A2                                                      
      WRITE(6,102)                                                      
 102  FORMAT(1X,' ELAB (MEV) =')                                        
      READ(5,*) EP              
      WRITE(6,103)
 103  FORMAT(1X,' TH-gamma (degree) =')
      READ(5,*) THGA
      WRITE(6,104)
 104  FORMAT(1X,' PHI-gamma (degree) =')
      READ(5,*) PHIGA
      WRITE(6,105)
 105  FORMAT(1X,' TH-lab-min (degree) =')
      READ(5,*) TH1
      WRITE(6,106)
 106  FORMAT(1X,' TH-lab-max (degree) =')
      READ(5,*) TH2
C
      VCM=0.04634*(A1/(A1+A2))*SQRT(EP/A1) 
C
      WRITE(6,108) A1,A2,EP,VCM,THGA,PHIGA                               
 108  FORMAT(1X,'A1= ',F4.0,' A2= ',F4.0,      
     1' EP = ',F7.2,'(MEV) V-CM/C = ',F10.5/
     21X,'TH-gamma = ',F7.2,'(degree) PHI-gamma = ',F7.2,'(degree)'/)         
C     
      WRITE(6,109)                                                      
 109  FORMAT(5X,'TH-lab    TH-CM    V-1/c',        
     1'     E-gamma   E-gamma   E-gamma   E-gamma'/
     224X,'V-2/c     0-deg     90-deg    180-deg   270-deg'/) 
C     
      DELTA=5.
      THGA=THGA*PI/180.
      PHIGA=PHIGA*PI/180.
C
   1  THETA=TH1*PI/180.
C
C
C
      PHI=0.*PI/180.
      CALL DOPP1
      EGAMMA11=EGAMMA
C
      PHI=90.*PI/180.
      CALL DOPP1
      EGAMMA21=EGAMMA
C
      PHI=180.*PI/180.
      CALL DOPP1
      EGAMMA31=EGAMMA
C
      PHI=270.*PI/180.
      CALL DOPP1
      EGAMMA41=EGAMMA
C      
      WRITE(6,110) TH1,THCM*180./PI,V1,EGAMMA11,EGAMMA21,EGAMMA31,EGAMMA41
 110  FORMAT(1X,2F10.3,5F10.5)
C
C
C
      PHI=0.*PI/180.
      CALL DOPP2
      EGAMMA12=EGAMMA
C
      PHI=90.*PI/180.
      CALL DOPP2
      EGAMMA22=EGAMMA
C
      PHI=180.*PI/180.
      CALL DOPP2
      EGAMMA32=EGAMMA
C
      PHI=270.*PI/180.
      CALL DOPP2
      EGAMMA42=EGAMMA
C      
      WRITE(6,110) TH1,THCM*180./PI,V2,EGAMMA12,EGAMMA22,EGAMMA32,EGAMMA42
C
C
C
      PHI=0.*PI/180.
      CALL DOPP3
      EGAMMA13=EGAMMA
C
      PHI=90.*PI/180.
      CALL DOPP3
      EGAMMA23=EGAMMA
C
      PHI=180.*PI/180.
      CALL DOPP3
      EGAMMA33=EGAMMA
C
      PHI=270.*PI/180.
      CALL DOPP3
      EGAMMA43=EGAMMA
C      
      WRITE(6,110) TH1,THCM*180./PI,V2,EGAMMA13,EGAMMA23,EGAMMA33,EGAMMA43
C
C
C
      PHI=0.*PI/180.
      CALL DOPP4
      EGAMMA14=EGAMMA
C
      PHI=90.*PI/180.
      CALL DOPP4
      EGAMMA24=EGAMMA
C
      PHI=180.*PI/180.
      CALL DOPP4
      EGAMMA34=EGAMMA
C
      PHI=270.*PI/180.
      CALL DOPP4
      EGAMMA44=EGAMMA
C      
      WRITE(6,110) TH1,THCM*180./PI,V1,EGAMMA14,EGAMMA24,EGAMMA34,EGAMMA44
C
C
C
      WRITE(21,209) TH1,
     11.454*EGAMMA11,1.454*EGAMMA21,1.454*EGAMMA31,1.454*EGAMMA41
      WRITE(22,209) TH1,
     11.257*EGAMMA12,1.257*EGAMMA22,1.257*EGAMMA32,1.257*EGAMMA42
      WRITE(23,209) TH1,
     11.257*EGAMMA13,1.257*EGAMMA23,1.257*EGAMMA33,1.257*EGAMMA43
      WRITE(24,209) TH1,
     11.454*EGAMMA14,1.454*EGAMMA24,1.454*EGAMMA34,1.454*EGAMMA44
 209  FORMAT(1X,5F10.5)                                        
C
      TH1=TH1+DELTA
      IF(TH1.GT.TH2) GOTO 3
      GOTO 1
C
    3 CONTINUE   
C
C
      write(21,221) THGA*180./PI 
      write(22,221) THGA*180./PI 
      write(23,221) THGA*180./PI 
      write(24,221) THGA*180./PI 
  221 format(10HW:th^ga^$:,F5.1,2H0@,
     132H/X(18) Y(1.15) s(6) POS(0.5,0.5))
      write(21,222) PHIGA*180./PI
      write(22,222) PHIGA*180./PI
      write(23,222) PHIGA*180./PI
      write(24,222) PHIGA*180./PI
  222 format(10HW:ph^ga^$:,F5.1,2H0@,
     132H/X(18) Y(1.12) s(6) POS(0.5,0.5))
C
C
C
      write(21,230)
      write(22,230)
      write(23,230)
      write(24,230)
  230 format(42HP: XMIN(10.) XMAX(50.) YMIN(1.1) YMAX(1.6)) 
      write(21,231)
      write(22,231)
      write(23,231)
      write(24,231)
  231 format(38HP: SHIF(0) SYMB(6) SCAL(0.6) FORM(1.2))
      write(21,232)
      write(22,232)
      write(23,232)
      write(24,232)
  232 format(35HP: POLY(2) AXFA(1.2) CMOD(3) UN(MM))
      write(21,233)
      write(22,233)
      write(23,233)
      write(24,233)
  233 format(34HP: XIN YIN XLI YLI NOSP TICI NODSN)
      write(21,234)
      write(22,234)
      write(23,234)
      write(24,234)
  234 format(9HP: WEI(0))
C
C
C                                                      
      STOP                                                             
      END                                                              
C 
      SUBROUTINE DOPP1
C                                                                       
      common/input/A1,A2,EP,THGA,PHIGA,VCM,THETA,PHI
      common/output/THCM,V1,V2,COSPG1,COSPG2,COSTG1,COSTG2,EGAMMA
C
C     projectile measured, projectile excited
C 
      PI=3.1415927                                                      
C
      THCM=THETA+ASIN((A1/A2)*SIN(THETA))
C             
      V1=VCM*SQRT(1+(A2/A1)**2+2*(A2/A1)*COS(THCM))
C
      COSPG1=COS(PHIGA)*COS(PHI)+SIN(PHIGA)*SIN(PHI)
C
      COSTG1=COS(THGA)*COS(THETA)+SIN(THGA)*SIN(THETA)*COSPG1
C
      EGAMMA=SQRT(1-V1**2)/(1.-V1*COSTG1)
C                                                      
      RETURN  
      END
C 
      SUBROUTINE DOPP2
C                                                                       
      common/input/A1,A2,EP,THGA,PHIGA,VCM,THETA,PHI
      common/output/THCM,V1,V2,COSPG1,COSPG2,COSTG1,COSTG2,EGAMMA
C
C     projectile measured, target excited
C 
      PI=3.1415927                                                      
C
      THCM=THETA+ASIN((A1/A2)*SIN(THETA))
C
      TH2=0.5*(PI-THCM)
C             
      V2=2.*VCM*COS(TH2)
C
      COSPG1=COS(PHIGA)*COS(PHI)+SIN(PHIGA)*SIN(PHI)
C
      COSTG2=COS(THGA)*COS(TH2)-SIN(THGA)*SIN(TH2)*COSPG1
C
      EGAMMA=SQRT(1-V2**2)/(1.-V2*COSTG2)
C                                                      
      RETURN     
      END
C 
      SUBROUTINE DOPP3
C
C     target nucleus measured, target excited
C                                                                       
      common/input/A1,A2,EP,THGA,PHIGA,VCM,THETA,PHI
      common/output/THCM,V1,V2,COSPG1,COSPG2,COSTG1,COSTG2,EGAMMA
C 
      PI=3.1415927                                  
C
      THCM=PI-2.*THETA
C
      V2=2.*VCM*COS(THETA)
C
      COSPG2=COS(PHIGA)*COS(PHI)+SIN(PHIGA)*SIN(PHI)
C
      COSTG2=COS(THGA)*COS(THETA)+SIN(THGA)*SIN(THETA)*COSPG2
C
      EGAMMA=SQRT(1-V2**2)/(1.-V2*COSTG2)
C                                                      
      RETURN                                                             
      END                                                              
C 
      SUBROUTINE DOPP4
C                                                                       
      common/input/A1,A2,EP,THGA,PHIGA,VCM,THETA,PHI
      common/output/THCM,V1,V2,COSPG1,COSPG2,COSTG1,COSTG2,EGAMMA
C
C     target nucleus measured, projectile excited
C 
      PI=3.1415927                                                      
C
      THCM=PI-2.*THETA
C             
      V1=VCM*SQRT(1+(A2/A1)**2+2.*(A2/A1)*COS(THCM))
C
      TH1=ACOS((VCM/V1)*(1.+(A2/A1)*COS(THCM)))
C
      COSPG2=COS(PHIGA)*COS(PHI)+SIN(PHIGA)*SIN(PHI)
C
      COSTG1=COS(THGA)*COS(TH1)-SIN(THGA)*SIN(TH1)*COSPG2
C
      EGAMMA=SQRT(1-V1**2)/(1.-V1*COSTG1)
C                                                      
      RETURN                                                             
      END                                                              
