program OpticalBench
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!This program simulates an optical system for frequency-dependent squeezing!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	use optocad
	use rsplot
!!!!Variable declaration, allocation, and initialization
	character(len=200),dimension(30,18) :: K = ''
	real(8) :: pst1(7) = (/2.,1.,1.e-1,1.e-2,1.e-3,1.e-4,1.e-5/),&
			   pst2(7) = (/2.e-1,1.e-1,1.e-2,1.e-3,1.e-4,1.e-5,1.e-6/),&
			   shift2, angle
	integer i, j
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!Character arrays describe optical components (all components in this program, no external files)
!!!!Basically the order of components is consistent with beam propagation
!!!!From Main Laser to BS1
    K(1,1)='b               .0, .0, w=.015, ag=0., p=2. ,z=22.          # Input beam'
	K(2,1)='d            t, 2.5, 0., 2., ag=0., m=1, r=0., t=1. # WPQ1  @cb0,-3'
    K(3,1)=' +           t, dm=0.1'
	K(4,1)='d            t, 5., 0., 2., ag=0., m=1, r=0., t=1. # WPH1  @cb0,3'
    K(5,1)=' +           t, dm=0.1'
    K(6,1)='d       	  t, 8., 0., 2., ag=0., m=1, r=0., t=1., f=10.0 # L1  @cb0,-3'
    K(7,1)=' +           t, dm=0.5'
	K(8,1)='d            t, 12., 0., 1., ag=0., m=1, r=0., t=1.'
    K(9,1)=' +           t, dm=1.'
	K(10,1)='d            t, 16.5, 0., 1., ag=0., m=5, r=0., t=1. # FI1  @cb2,-3'
    K(11,1)=' +           t, dm=3.'
	K(12,1)='d            t, 23., 0., 1., ag=0., m=1, r=0., t=1.'
    K(13,1)=' +           t, dm=1.'
	K(14,1)='d            t, 25., 0., 2., ag=0., m=1, r=0., t=1. # WPH2  @cb0,3'
    K(15,1)=' +           t, dm=0.1'
	K(16,1)='d            t, 27.5, 0., 1., ag=0., m=5, r=0., t=1. # EOM(SHG) @cb2,-3'
    K(17,1)=' +           t, dm=3.'
	K(18,1)='d       	srt, 42.5, 0., 2., ag=-45., m=1, r=.9, t=.1 # BS1  @cb1,-1'
    K(19,1)=' +           t, dm=0.5'
!!!!Reflection of BS1 to SHG input
    K(1,2)='d            t, 42.5, 29., 2., ag=90., m=1, r=0., t=1., f=5.02 #L2 @cb3,0'
    K(2,2)=' +           t, dm=0.5'
    K(3,2)='d            r, 42.5, 35., 2., ag=45., m=1, r=1., t=0. #M2 @cb-3,1'
    K(4,2)=' +           t, dm=0.5'
    K(5,2)='d       	  t, 37.5, 35., 2., ag=180., m=1, r=0., t=1., f=6.29 #L3  @cb-1,-3'
    K(6,2)=' +           t, dm=0.5'
	K(7,2)='d            t, 32.5, 35., 2., ag=180., m=1, r=0., t=1. # WPH3 @cb0,3'
    K(8,2)=' +           t, dm=0.1'
    K(9,2)='d 	          t, 24.5, 35., 2., ag=180., m=1, r=0., t=1., f=12.5 # L4   @cb-1,-3'
    K(10,2)=' +           t, dm=0.5'
    K(11,2)='d           rh, 20., 35., 2., ag=135., m=1, r=1., t=0. # M3  @cb-1,-1'
    K(12,2)=' +           t, dm=0.5'
	K(13,2)='d      strstrh, 20., 25., 2., ag=-135., m=1, r=.99, t=.01 # DM4  @cb0,-3'
    K(14,2)=' +           t, dm=0.5'
	K(15,2)='c            d, 20., 40., 2., ag=90. #Damper @0,2' !Beam damper
!!!!SHG cavity and its transmission
	K(1,3)='d      crtstrh,  24.5, 25., 2., ag=180., m=1, r=0.92, t=0.08, rc=-2.5 # SHG @cb0,3'
    K(2,3)=' +           t, dm=.5'
	K(3,3)='d            t,  27.2, 25., 1., ag=.0, m=4, r=.0, t=1. # LN @cb0,-3'
    K(4,3)=' +         str, dm=0.65, r=0.999, t=0.001, rc=1.2'
	K(5,3)='d            r, 32.5, 25., 2., ag=45., m=1, r=1., t=0. # M5  @cb3,2'
    K(6,3)=' +           t, dm=0.5'
	K(7,3)='d            t, 32.5, 20., 2., ag=-90., m=1, r=0., t=1., f=10. # L5  @cb3,0'
    K(8,3)=' +           t, dm=0.5'
	K(9,3)='d            r, 32.5, 15., 2., ag=-45., m=1, r=1., t=0. # M6  @cb2,-1'
    K(10,3)=' +           t, dm=0.5'
	K(11,3)='c            d, 15., 15., 2., ag=0.'
	K(12,3)='c            d, 20., 20., 2., ag=-90. #Damper @0,-1' !Beam damper
!!!!Green beam starting point (in Lithium Niobate)
    K(1,4)='x                lam=532e-7'
    K(2,4)='b               27.5, 25., .00265, ag=180.,p=0.15             # Input beam'
!!!!SHG output (green) to Mach-Zehnder interferometer
!!!!Simulating a dichroic mirror is triky (be careful to difference of str and srt!)
	K(1,5)='d            t, 16., 25., 2., ag=180., m=1, r=0., t=1., f=20. # L6  @cb0,-3'
    K(2,5)=' +           t, dm=0.5'
	K(3,5)='c            d, 11.5, 20., 2., ag=-90. #Damper @0,-1' !Beam damper
	K(4,5)='c            d, 7.5, 30., 2., ag=90. #Damper @0,-1' !Beam damper
	K(5,5)='d          srt, 11.25, 25., 2., ag=135., m=1, r=.01, t=.99 # DM7  @cb-1,2'
    K(6,5)=' +           t, dm=0.5'
	K(7,5)='c         strh, 11.25, 25., 2., ag=135., m=1, r=.99, t=.01 # DM7ir'
	K(8,5)='d          srt, 7.5, 25., 2., ag=-135., m=1, r=.01, t=.99 # DM8  @cb-1,-2'
    K(9,5)=' +           t, dm=0.5'
	K(10,5)='c           rh, 7.5, 25., 2., ag=-135., m=1, r=1., t=.0 # DM8ir'
	K(11,5)='d            t, 0., 25., 1., ag=180., m=1, r=0., t=1.'
    K(12,5)=' +           t, dm=1.'
	K(13,5)='d            t, -2., 25., 1., ag=180., m=5, r=0., t=1. # FI2 @cb0,-3'
    K(14,5)=' +           t, dm=2.'
	K(15,5)='d            t, -5., 25., 1., ag=180., m=1, r=0., t=1.'
    K(16,5)=' +           t, dm=1.'
	K(17,5)='d            t, -13., 25., 2., ag=0., m=1, r=0., t=1. # WPH4 @cb0,3'
    K(18,5)=' +           t, dm=0.1'
	K(19,5)='d            r, -17.5, 25., 2., ag=-135., m=1, r=1., t=0. # M9  @cb-1,-1'
    K(20,5)=' +           t, dm=0.5'
	K(21,5)='d            t, -17.5, 41., 2., ag=90., m=1, r=0., t=1., f=20. # L7  @cb-1,-1'
    K(22,5)=' +           t, dm=0.5'
	K(23,5)='d            t, -17.5, 55., 1., ag=90., m=5, r=0., t=1. # EOM(FC) @cb4,0'
    K(24,5)=' +           t, dm=3.'
!!!!Mach-Zehnder (MZ)
	K(1,6)='d          str, -17.5, 70., 2., ag=135., r=0.5, t=0.5 #MZ @-1,-1' !Input BS
	K(2,6)=' +           t, dm=0.5'
	K(3,6)='d       	  r, -17.5-shift2, 75.+shift2, 2., ag=135., r=1., t=0.' !Mirror for transmitted beam
	K(4,6)=' +           t, dm=0.5'
	K(5,6)='d            r, -12.5, 70., 2., ag=-45., r=1., t=0.' !Mirror for reflected beam
	K(6,6)=' +           t, dm=0.5'
	K(7,6)='d     i(rt)srt, -12.5, 75., 2., ag=135., r=0.5, t=0.5' !Output BS
	K(8,6)=' +           t, dm=0.5'
	K(9,6)='c            d, -7.5, 75., 2., ag=0. #PD @0,-1' !Beam damper
!!!!After MZ to Mode Cleaner (Green)
	K(1,7)='d            t, -12.5-shift, 77.5, 2., ag=90., m=1, r=0., t=1., f=10. #L8 @3,0'
	K(2,7)=' +           t, dm=0.5'
	K(3,7)='d          str, -12.5, 80.+shift2, 2., ag=135., m=1, r=.5, t=.5 # BS @cb-3,0'
    K(4,7)=' +           t, dm=.5'
	K(5,7)='d            t,  .0, 80., 2., ag=0.0, m=1, r=0., t=1., f=10.0 # L11  @cb0,3'
    K(6,7)=' +           t, dm=0.5'	
	K(7,7)='d            t, 13.5, 80., 2., ag=0.0, m=1, r=0., t=1., f=5. # L12  @cb0,3'
    K(8,7)=' +           t, dm=0.5'
	K(9,7)='d       	  r, 17.5, 80., 2., ag=45., r=1., t=0. # M11 @cb3,3'
	K(10,7)=' +           t, dm=0.5'
	K(11,7)='d       	  r, 17.5, 57.5, 2., ag=-135., r=1., t=0. # M12 @cb-3,-3' 
	K(12,7)=' +           t, dm=0.5'
!!!!Mode Cleaner (Green) Cavity
	K(1,8)='d       crtstr, 25., 57.5+shift, 1.5, ag=-135.-angle+1.1, m=1, r=0.98, t=0.02 # MC(Green) @cb2,-4'
    K(2,8)=' +           t, dm=0.5'
	K(3,8)='d          str, 27.2, 57.5+shift, 1.5, ag=-45.+angle-1.1, m=1, r=0.98, t=0.02 '
    K(4,8)=' +           t, dm=0.5'
	K(5,8)='d            r, 26.1, 85.6+shift, 1.5, ag=90., m=1, r=.99, t=.01, rc=-100.'
    K(6,8)=' +           t, dm=0.5'
!!!!After Mode Cleaner (Green) to OPO
	K(1,9)='d       	srt, 40., 57.5, 2., ag=45., r=.9, t=.1 # BS9 @cb3,3'
	K(2,9)=' +           t, dm=0.5'
	K(3,9)='d       	  r, 47.5, 57.5, 2., ag=-45., r=1., t=0. # M23 @cb-3,-3' 
	K(4,9)=' +           t, dm=0.5'
	K(5,9)='d            t, 47.5, 62.5, 2., ag=90., m=1, r=0., t=1., f=5. # L14  @cb0,3'
    K(6,9)=' +           t, dm=0.5'
	K(7,9)='c            d, 47.5, 67.5, 2., ag=90. #PD @cb2,0'
	K(8,9)='d       	  r, 40., 50., 2., ag=-135., r=1., t=0. # M24 @cb-3,-3' 
	K(9,9)=' +           t, dm=0.5'
	K(10,9)='d       	  r, 95., 50., 2., ag=-45., r=1., t=0. # M25 @cb-3,-3' 
	K(11,9)=' +           t, dm=0.5'
	K(12,9)='d       	  r, 24.5, 50., 2., ag=-46., r=1., t=0. # M26 @cb0,-3' 
	K(13,9)=' +           t, dm=0.5'
	K(14,9)='d            t, 15., 50., 2., ag=-180., m=1, r=0., t=1., f=5. # L15  @cb0,-3'
    K(15,9)=' +           t, dm=0.5'
	K(16,9)='c            d, 10., 50., 2., ag=-180. #PD @cb0,-2'
!!!!After MZ to vacuume chamber
	K(1,10)='d            r, -12.5-2*shift2, 97.5, 2., ag=135., m=1, r=1., t=0. # M10 @cb-3,0'
    K(2,10)=' +           t, dm=.5'
	K(3,10)='d            t, -5., 97.5, 2., ag=0., m=1, r=0., t=1., f=10.# L9 @cb0,-3'
    K(4,10)=' +           t, dm=.5'
	K(5,10)='d            t, 25., 97.5, 2., ag=0., m=1, r=0., t=1., f=20.# L10 @cb0,-3'
    K(6,10)=' +           t, dm=.5'
	K(7,10)='d            t, 32.5, 97.5, 1., ag=0., m=1, r=0., t=1.'
    K(8,10)=' +           t, dm=1.'
	K(9,10)='d            t, 34.5, 97.5, 1., ag=0., m=5, r=0., t=1. # FI3 @cb0,-3'
    K(10,10)=' +           t, dm=2.'
	K(11,10)='d            t, 37.5, 97.5, 1., ag=0., m=1, r=0., t=1.'
    K(12,10)=' +           t, dm=1.'
	K(13,10)='d            t, 40., 97.5, 2., ag=0., m=1, r=0., t=1. # WPH5 @cb0,3'
    K(14,10)=' +           t, dm=0.1'
	K(15,10)='d       	   r, 60., 97.5, 2., ag=-45., r=1., t=0. # M13 @cb0,-3'
	K(16,10)=' +           t, dm=0.5'
	K(17,10)='d       	   t, -12.5-2*shift, 91.75, 1.5, ag=90., m=1, r=0., t=1. #AOMblock'
	K(18,10)=' +           t, dm=1.5'
	K(19,10)='c      {st-1t}, -12.5-2*shift, 92.5, 1.5, ag=90., t=.1, t-1=.9, gf=200. #AOM @cb3,0'
!!!!Transmission of BS1
	K(1,11)='d       	   t, 75.5, 0.-shift, 2., ag=0., m=1, r=0., t=1., f=50. # L13 @cb3,0'
	K(2,11)=' +           t, dm=0.5'
	K(3,11)='d       	 str, 102.5, 0., 2., ag=-45., m=1, r=.1, t=.9 # BS3  @cb1,-3'
    K(4,11)=' +           t, dm=0.5'
!!!!Reflection of BS3
	K(5,11)='d       	 str, 102.5, 10., 2., ag=135., m=1, r=.5, t=.5 # BS4  @cb1,-1'
    K(6,11)=' +           t, dm=0.5'
!!!!Transmission of BS4
	K(7,11)='d       	   r, 102.5, 25.+2*shift2, 2., ag=135., r=1., t=0. # M17 @cb0,-3'
	K(8,11)=' +           t, dm=0.5'
	K(9,11)='d       	 str, 107.5, 25., 2., ag=45., m=1, r=.1, t=.9 # BS10  @cb1,-1'
    K(10,11)=' +           t, dm=0.5'
	K(11,11)='c            d, 107.5, 20., 2., ag=-90. #PD @cb2,0'
	K(12,11)='c            d, 112.5, 25., 2., ag=0. #PLL @cb2,0'
!!!!Transmission of BS3
	K(13,11)='d       	   r, 142.5, .0, 2., ag=-45., r=1., t=0. # M14 @cb0,-3'
	K(14,11)=' +           t, dm=0.5'
	K(15,11)='d       	   t, 142.5, 30., 2., ag=90., r=0., t=1., f=100. # L14 @cb3,0'
	K(16,11)=' +           t, dm=0.5'
	K(17,11)='d       	   t, 142.5, 60., 2., ag=90., r=0., t=1., f=100. # L15 @cb3,0'
	K(18,11)=' +           t, dm=0.5'
	K(19,11)='d       	   r, 142.5, 95., 2., ag=45., r=1., t=0. # M15 @cb0,-3'
	K(20,11)=' +           t, dm=0.5'
	K(21,11)='d       	   r, 65., 95., 2., ag=-135., r=1., t=0. # M16 @cb0,-3'
	K(22,11)=' +           t, dm=0.5'
!!!!Aux. Laser 2 beam start
    K(1,12)='o            xo=75., yo=15.' !Set an offset
    K(2,12)='b            .0, 0., w=.018, ag=0., z=20.                       # Input beam'
!!!!Aux2 beam path
	K(3,12)='d            t, 2.5, 0., 2., ag=0., m=1, r=0., t=1. # WPQ  @cb0,-3'
    K(4,12)=' +           t, dm=0.1'
	K(5,12)='d            t, 5., 0., 2., ag=0., m=1, r=0., t=1. # WPH  @cb0,3'
    K(6,12)=' +           t, dm=0.1'
    K(7,12)='d       	 t, 10., 0., 2., ag=0., m=1, r=0., t=1., f=30.0 # L  @cb0,-3'
    K(8,12)=' +           t, dm=0.5'
	K(9,12)='d            t, 12., 0., 1., ag=0., m=1, r=0., t=1.'
    K(10,12)=' +           t, dm=1.'
	K(11,12)='d            t, 16.5, 0., 1., ag=0., m=5, r=0., t=1. # FI  @cb2,-3'
    K(12,12)=' +           t, dm=3.'
	K(13,12)='d            t, 23., 0., 1., ag=0., m=1, r=0., t=1.'
    K(14,12)=' +           t, dm=1.'
	K(15,12)='d            t, 25., 0., 2., ag=0., m=1, r=0., t=1. # WPH  @cb0,3'
    K(16,12)=' +           t, dm=0.1'
	K(17,12)='d       	 str, 32.5, 0., 2., ag=45., m=1, r=.5, t=.5 # BS5  @cb1,-1'
    K(18,12)=' +           t, dm=0.5'
!!!!Reflection from BS5
	K(19,12)='d       	 str, 32.5, -5., 2., ag=-135., m=1, r=.5, t=.5 # BS6  @cb1,-1'
    K(20,12)=' +           t, dm=0.5'
!!!!Reflection from BS6 to PLL
	K(21,12)='c            d, 37.5, -5., 2., ag=0. #PLL @0,-1' !PLL
!!!!Transmission from BS6 to PD
	K(22,12)='c            d, 32.5, -10., 2., ag=-90. #PD @0,-1' !PD
!!!!Transmission from BS5 to OPO
	K(23,12)='d       	   r, 50., 0., 2., ag=-45., r=1., t=0. # M17 @cb0,-3'
	K(24,12)=' +           t, dm=0.5'
	K(25,12)='d            t, 50., 10., 1., ag=90., m=5, r=0., t=1. # EOM(OPO) @cb4,0'
    K(26,12)=' +           t, dm=3.'
!!!!Polarizing Beam Splitter 1
    K(1,13)='o            xo=75., yo=15.' !Set an offset
	K(2,13)='c       	  t, 50., 22.5-2./1.414, rd=2./1.414, ag=90., m=1, r=.0, t=1. # PBS  @cb3,-2' !Ploting PBS cube
	K(3,13)=' +      	  t, rd=2./1.414, ag=180., m=1, r=.0, t=1.'
	K(4,13)=' +      	  t, rd=2./1.414, ag=270., m=1, r=.0, t=1.'
	K(5,13)=' +       	  t, rd=2./1.414, ag=0., m=1, r=.0, t=1.'
	K(6,13)='c       	srt, 50., 22.5, 2., ag=-45., r=.995, t=.005 #s-pol' !Surface for s-pol.(Aux1)
	K(7,13)='c        srth, 50., 22.5, 2., ag=-45., r=.005, t=.995 #p-pol' !Surface for p-pol.(Aux2)
!!After PBS1 to OPO
    K(1,14)='o            xo=75., yo=15.' !Set an offset
	K(2,14)='d       	  r, 50., 35., 2., ag=45., r=1., t=0. # M18 @cb-3,-3' 
	K(3,14)=' +           t, dm=0.5'
	K(4,14)='d       	  r, 30., 35., 2., ag=-135., r=1., t=0. # M19 @cb-3,-3' 
	K(5,14)=' +           t, dm=0.5'
	K(6,14)='d       	srtr, 30., 52.5, 2., ag=45., r=.1, t=.9 # BS7 @cb2,2' 
	K(7,14)=' +           t, dm=0.5'
	K(8,14)='c            d, 25., 52.5, 2., ag=180. #BD @-1,1' !Beam damper
	K(9,14)='d       	  r, 30., 57.5+2*shift2, 2., ag=45., r=1., t=0. # M20 @cb2,2' 
	K(10,14)=' +           t, dm=0.5'
	K(11,14)='d       	  r, 20., 57.5, 2., ag=135., r=1., t=0. # M21 @cb-2,2' 
	K(12,14)=' +           t, dm=0.5'
	K(13,14)='c            d, 35., 42.5, 2., ag=-90. #PD(OPO Phase) @0,-1'
!!!!Polarizing Beam Splitter 2
	K(1,15)='c       	  t, 110., 67.5-2./1.414, rd=2./1.414, ag=90., m=1, r=.0, t=1. # PBS  @cb2,2' !Ploting PBS cube
	K(2,15)=' +      	  t, rd=2./1.414, ag=180., m=1, r=.0, t=1.'
	K(3,15)=' +      	  t, rd=2./1.414, ag=270., m=1, r=.0, t=1.'
	K(4,15)=' +       	  t, rd=2./1.414, ag=0., m=1, r=.0, t=1.'
	K(5,15)='c       	srt, 110., 67.5, 2., ag=45., r=.995, t=.005 #s-pol' !Surface for s-pol.(Aux1)
	K(6,15)='c        srth, 110., 67.5, 2., ag=45., r=.005, t=.995 #p-pol' !Surface for p-pol.(Aux2)
!!!!Aux. Laser 1 beam start
    K(1,16)='o            xo=75., yo=37.5' !Set an offset
    K(2,16)='b            0., 0., w=.019, ag=0., z=18.                       # Input beam'
!!!!Aux1 beam path
	K(3,16)='d            t, 2.5, 0., 2., ag=0., m=1, r=0., t=1. # WPQ  @cb0,-3'
    K(4,16)=' +           t, dm=0.1'
	K(5,16)='d            t, 5., 0., 2., ag=0., m=1, r=0., t=1. # WPH  @cb0,3'
    K(6,16)=' +           t, dm=0.1'
    K(7,16)='d       	 t, 8., 0., 2., ag=0., m=1, r=0., t=1., f=30.0 # L  @cb0,-3'
    K(8,16)=' +           t, dm=0.5'
	K(9,16)='d            t, 12., 0., 1., ag=0., m=1, r=0., t=1.'
    K(10,16)=' +           t, dm=1.'
	K(11,16)='d            t, 16.5, 0., 1., ag=0., m=5, r=0., t=1. # FI  @cb2,-3'
    K(12,16)=' +           t, dm=3.'
	K(13,16)='d            t, 23., 0., 1., ag=0., m=1, r=0., t=1.'
    K(14,16)=' +           t, dm=1.'
	K(15,16)='d            t, 25., 0., 2., ag=0., m=1, r=0., t=1. # WPH  @cb0,3'
    K(16,16)=' +           t, dm=0.1'
	K(17,16)='d       	 str, 32.5+shift2, 0., 2., ag=45., m=1, r=.5, t=.5 # BS8  @cb1,-1'
    K(18,16)=' +           t, dm=0.5'
!!!!OPO Cavity
	K(1,17)='d       crtstr,  95., 67.5, 2., ag=90., m=1, r=.92, t=.08, rc=-2.5 # OPO @cb-3,0'
    K(2,17)=' +           t, dm=.5'
	K(3,17)='c           rh,  95., 67.5, 2., ag=90., m=1, r=.92, t=.08, rc=-2.5 #green'
	K(4,17)='d            t,  95., 64.8, 1., ag=.-90, m=4, r=.0, t=1. # LN @cb0,-3'
    K(5,17)=' +      crtstr, dm=0.65, r=0.999, t=0.001, rc=1.2'
!!!!OPO Output
	K(1,18)='d             hr, 95., 57.5, 2., ag=-45., r=1., t=0. # DM1 @cb3,0' 
	K(2,18)=' +            h, dm=0.5'
	K(3,18)='c             t, 95., 57.5, 2., ag=135., r=0., t=1. # DM1g' 
	K(4,18)='d       	  h, 87.5, 57.5, 2., ag=-45., r=1., t=0. # FM @cb0,-3' 
	K(5,18)=' +            h, dm=0.5'
	K(6,18)='d       	  h, 82.5, 57.5, 2., ag=-135., r=1., t=0. # FM @cb0,-3' 
	K(7,18)=' +            h, dm=0.5'
	K(8,18)='d       	  r, 65., 57.5, 2., ag=-135., r=1., t=0. #M27' 
	K(9,18)=' +            t, dm=0.5'
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!OptoCad initialization and setups
    call oc_init(unit='cm',pform='A0_L')
    call oc_frame(-26.3,-13.7,153.7,106.3,50.,70.,.6,glp=0)
    call oc_set(n=1, fslb=6., &
				rix=(/1.4496_8, 1.4285_8, 1.5066_8, 2.2321, 1.954_8/), & 
				! Fused Silica, CaF2, N-BK7, LiNb, TGG? (https://refractiveindex.info/)
				print='rs s2 act w0t z1t ang Ang rd dx dy pw lb')
!!!!PS PLOT drawing section
	call ps_quad(frm=1,xu=-26.3,yu=-13.7,wu=180.,hu=5.05,fill=14) !Post bottom
	call ps_quad(frm=1,xu=153.7 - 5.05,yu=-13.7,wu=5.05,hu=120.,fill=14) !Post right
	call ps_quad(frm=1,xu=-26.3,yu=-13.7,wu=5.05,hu=120.,fill=14) !Post left
	call ps_quad(frm=1,xu=-26.3,yu=106.3 - 5.05,wu=180.,hu=5.05,fill=14) !Post top
	call ps_quad(frm=1,xu=-19.2,yu=-7.45,wu=19.2,hu=14.9,fill=13) !Master laser
	call ps_text(frm=1,str='Main Laser',xu=-19.2,yu=8.,fs=6.) 
	call ps_quad(frm=1,xu=75.-19.2,yu=15.-14.9/2,wu=19.2,hu=14.9,fill=13) !Aux laser 1
	call ps_text(frm=1,str='Aux. Laser 2',xu=55.,yu=23.,fs=6.) 
	call ps_quad(frm=1,xu=75.-19.2,yu=37.5-14.9/2,wu=19.2,hu=14.9,fill=13) !Aux laser 2
	call ps_text(frm=1,str='Aux. Laser 1',xu=55.,yu=46.,fs=6.)
	call ps_quad(frm=1,xu=12.,yu=-2.,wu=12.,hu=4.,fill=14) !FI1
	call ps_quad(frm=1,xu=27.,yu=-2.,wu=4.,hu=4.,fill=14) !EOM(SHG)
	call ps_quad(frm=1,xu=23.5,yu=21.5,wu=5.,hu=7.,fill=14) !SHG housing
	call ps_quad(frm=1,xu=5.,yu=12.5,wu=10.,hu=5.,fill=14) !PD(SHG)
	call ps_text(frm=1,str='PD(SHG)',xu=5.,yu=10.,fs=6.) 
	call ps_quad(frm=1,xu=-6.,yu=23.,wu=6.,hu=4.,fill=14) !FI2
	call ps_quad(frm=1,xu=-19.5,yu=54.5,wu=4.,hu=4.,fill=14) !EOM(FC)
	call ps_quad(frm=1,xu=123.,yu=24.5,wu=4.,hu=4.,fill=14) !EOM(OPO)
	call ps_quad(frm=1,xu=-14.5-2*(shift(45._8,1.4496_8,0.5_8) + 0.0125269),yu=91.5,wu=4.,hu=2.,fill=14) !AOM
	call ps_quad(frm=1,xu=21.5,yu=55.5,wu=9.0,hu=31.5,fill=14) !Mode Cleaner(Green)
	call ps_quad(frm=1,xu=32.5,yu=95.5,wu=6.,hu=4.,fill=14) !FI3
	call ps_quad(frm=1,xu=75.+12.,yu=15.-2.,wu=12.,hu=4.,fill=14) !FI Aux2
	call ps_quad(frm=1,xu=75.+12.,yu=37.5-2.,wu=12.,hu=4.,fill=14) !FI Aux1
	call ps_quad(frm=1,xu=95.-3.5,yu=63.5,wu=7.,hu=5.,fill=14) !OPO housing
	call ps_text(frm=1,str='Optical Bench',xl=0.,yt=10.,fs=25.,fn='H') 
!!!!"Shift" is a value used to compensate beam path shift caused by any optics.
	call oc_bind('shift',shift(45._8,1.4496_8,0.5_8))
	call oc_bind('shift2',shift(45._8,1.4496_8,0.5_8)*tan(45._8*3.14159265358979_8/180._8))
	call oc_bind('angle',atan(1.1_8/28.1_8)*180._8/3.14159265358979_8)
!!!!Reading optical components data from character array K and plotting surfaces.
	do i = 1, 18
		call oc_input(K(:,i))
	end do
	call oc_surf(lw=.1)
!!!!Beam color for 1.IR and 2.green.
	ico1=ps_color((/1.,.8,.8/))
	icm1=ps_color((/1.,.6,.6/))
	ici1=ps_color((/1.,.0,.0/))
	ico2=ps_color((/.8,1.,.8/))
	icm2=ps_color((/.6,1.,.6/))
	ici2=ps_color((/.0,1.,.0/))
!!!!Trace ray segments (order : Main => Green => Aux2 => Aux1)
!!!!Changing the plotting order may cause problems.
	do j = 1, 4	
		if (j /= 2) then
			call oc_trace(j)                      
			call oc_beam(4.,fill=ico1,pst=pst1)
			call oc_beam(2.4,fill=icm1,pst=pst1)
			call oc_beam(1.2,fill=ici1,pst=pst1)
			call oc_beam(.0,1,5,.1,rnc=1)
		else
			call oc_trace(j)                      
			call oc_beam(4.,fill=ico2,pst=pst2)
			call oc_beam(2.4,fill=icm2,pst=pst2)
			call oc_beam(1.2,fill=ici2,pst=pst2)
			call oc_beam(.0,1,5,.1,rnc=1)
		end if
	end do
!!!!Exit OptoCAD (Good job!)
    call oc_exit
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
contains
	real(8) function shift(theta1, n, a)
		!!!This function calculates beam path shift 
		!!!caused by optical components (mainly beam splitter) of refractive index n.
		!!!theta1[deg]: angle of incidence (normally 45 deg.).
		!!!theta2[rad]: angle of refraction (normally less than theta1).
		!!!n: refractive index of a component.
		!!!a[cm]: thickness of a component.
		real(8), intent(in) :: theta1, n, a
		real(8) pi, theta2
		pi = 3.14159265358979_8
		theta2 = asin(sin(theta1*pi/180._8)/n)
		shift = a*sin(theta1*pi/180._8) - a*cos(theta1*pi/180._8)*tan(theta2)
	end function shift
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
end program OpticalBench