param = {ay1 -> 20.27, ax2 -> 54.27,  
        rp -> 1.5 10^(-18),  lw -> 1.2,  
        Iw -> 120,  vl -> 3 10^8,  el -> 1.6 10^-19,  
        gamma0 ->  27.7}

coef  =  rp Iw lw/(2 Pi gamma0 el vl) //. param

bx1ave  =  (45.6 + 49.69)/2. ;

by1ave  =  (52.97 + 48.73)/2. ;

bx2ave  =  (51.91 + 56.36)/2.

by2ave  =  (46.61 + 42.76)/2.

Solve[(3 x^6 + 9 u^5+20) == 0, x]

Clear[ax1]

sol  =  
      Solve[((ay1 + dy) 2 - 
                dx^2) + (dy^2 - (ax2 - dx)^2) == 0, 
        dy]

sol  =  
        Solve[((ay1 + dy) 2 - 
                    dx^2) ((ax2 - dx)^2 + 
                      dy^2)^2 + (dy^2 - (ax2 - 
                          dx)^2) (dx^2 + (ay1 + dy)^2)^2 == 
            0, dy];

dy[dx_] := 
      dy //. FindRoot[((ay1 + dy)^2 - 
                      dx^2) ((ax2 - dx)^2 + 
                        dy^2)^2 + (dy^2 - (ax2 - 
                            dx)^2) (dx^2 + (ay1 + 
                              dy)^2)^2 == 0 //. param, {dy, 
            20}]

dy[16. ]

Qx[dy[15], 15] //. param

Qx[dyh_, dxh_] := 
      coef (bx1ave ((ay1 + dy[dxh])^2 - 
                    dxh^2)/((ay1 + dy[dxh])^2 + dxh^2)^2 + 
            bx2ave ((dy[dxh])^2 - (ax2 - 
                          dxh)^2)/((dy[dxh])^2 + (ax2 - 
                            dxh)^2)^2) 10^6

Qy[dyh_, dxh_] := 
      coef (-by1ave ((ay1 + dy[dxh])^2 - 
                    dxh^2)/((ay1 + dy[dxh])^2 + dxh^2)^2 - 
            by2ave ((dy[dxh])^2 - (ax2 - 
                          dxh)^2)/((dy[dxh])^2 + (ax2 - 
                            dxh)^2)^2) 10^6

kappa[dyh_, dxh_] := 
      coef 2 Abs[(Sqrt[
                    bx1ave by1ave] (2 (ay1 + 
                            dy[dxh]) dxh )/((ay1 + dy[dxh])^2 + 
                          dxh^2)^2 + 
                Sqrt[bx2ave by2ave] (2 dy[
                          dxh] (ax2 - 
                            dxh))/((dy[dxh])^2 + (ax2 - 
                                dxh)^2)^2) 10^6]

p1 = 
      Plot[Qx[dy[dx], dx] //. param, {dx, 5, 30}, 
        PlotRange -> {{5, 30}, {- .002, 0.002}}]

p2  =  
      Plot[Qy[dy[dx], dx] //. param, {dx, 5, 30}, 
        PlotRange -> {{5, 30}, {-0.002, 0.002}}]

Show[p1, p2]

p1kappa = 
      Plot[kappa[dy[dx], dx] //. param, {dx, 5, 30}, 
        PlotRange -> {{5, 30}, {0, 0.015}}]

Qx2[dyh_, dxh_] := 
      2.  coef (bx1ave ((ay1 + dy[dxh])^2 - 
                    dxh^2)/((ay1 + dy[dxh])^2 + dxh^2)^2 + 
            0.  ((dy[dxh])^2 - (ax2 - 
                          dxh)^2)/((dy[dxh])^2 + (ax2 - 
                            dxh)^2)^2) 10^6

Qy2[dyh_, dxh_] := 
      2.  coef (-by1ave ((ay1 + dy[dxh])^2 - 
                    dxh^2)/((ay1 + dy[dxh])^2 + dxh^2)^2 - 
            0.  ((dy[dxh])^2 - (ax2 - 
                          dxh)^2)/((dy[dxh])^2 + (ax2 - 
                            dxh)^2)^2) 10^6

Qx3[dyh_, dxh_] := 
      2.  coef (bx1ave ((ay1 + dy[dxh])^2 - 
                    dxh^2)/((ay1 + dy[dxh])^2 + dxh^2)^2 0.  + 
            bx2ave ((dy[dxh])^2 - (ax2 - 
                          dxh)^2)/((dy[dxh])^2 + (ax2 - 
                            dxh)^2)^2) 10^6

Qy3[dyh_, dxh_] := 
      2.  coef (-by1ave ((ay1 + dy[dxh])^2 - 
                    dxh^2)/((ay1 + dy[dxh])^2 + dxh^2)^2 0.  -  
            by2ave ((dy[dxh])^2 - (ax2 - 
                          dxh)^2)/((dy[dxh])^2 + (ax2 - 
                            dxh)^2)^2) 10^6

kappayy[dyh_, dxh_] := 
      coef 4  Abs[(Sqrt[
                    bx1ave by1ave] (2 (ay1 + 
                            dy[dxh]) dxh )/((ay1 + dy[dxh])^2 + 
                          dxh^2)^2 +  
                0 Sqrt[
                    bx2ave by2ave] (2 dy[
                          dxh] (ax2 - 
                            dxh))/((dy[dxh])^2 + (ax2 - 
                                dxh)^2)^2) 10^6]

kappaxx[dyh_, dxh_] := 
      coef 4  Abs[(Sqrt[
                    bx1ave by1ave] (2 (ay1 + 
                            dy[dxh]) dxh )/((ay1 + dy[dxh])^2 + 
                          dxh^2)^2 0 +   Sqrt[
                    bx2ave by2ave] (2 dy[
                          dxh] (ax2 - 
                            dxh))/((dy[dxh])^2 + (ax2 - 
                                dxh)^2)^2) 10^6]

Qx3[dy[10. ], 10. ] //. param

p3 = 
      Plot[Qx2[dy[dx], dx] //. param, {dx, 5, 30}, 
        PlotRange -> {{5, 30}, {-0.002, 0.002}},  
        PlotStyle -> Dashing[{0.02, 0.03}]]

p4  =  
      Plot[Qy2[dy[dx], dx] //. param, {dx, 5, 30}, 
        PlotRange -> {{5, 30}, {-0.002, 0.002}},  
        PlotStyle -> Dashing[{0.02, 0.03}]]

Show[p3, p4]

Show[p3, p1, p2, p4, 
      AxesLabel -> {"Dx [mm]", "Delta Q_x,y"}]

p3b = 
      Plot[Qx3[dy[dx], dx] //. param, {dx, 5, 30}, 
        PlotRange -> {{5, 30}, {-0.002, 0.002}},  
        PlotStyle -> Dashing[{0.02, 0.03}]]

p4b  =  
      Plot[Qy3[dy[dx], dx] //. param, {dx, 5, 30}, 
        PlotRange -> {{5, 30}, {-0.002, 0.002}},  
        PlotStyle -> Dashing[{0.02, 0.03}]]

Show[p3, p4, p3b, p4b]

Export["DetuningWithDx.eps", %]

Show[p3, p1, p2, p4, 
      AxesLabel -> {"Dx [mm]", "Delta Q_x,y"}]

p3kappa = 
      Plot[kappaxx[dy[dx], dx] //. param, {dx, 5, 30}, 
        PlotRange -> {{5, 30}, {-0.02, 0.02}},  
        PlotStyle -> Dashing[{0.02, 0.06}]]

p4kappa = 
      Plot[kappayy[dy[dx], dx] //. param, {dx, 5, 30}, 
        PlotRange -> {{5, 30}, {-0.02, 0.02}},  
        PlotStyle -> Dashing[{0.02, 0.03}]]

Show[p1kappa, p3kappa, p4kappa, 
      AxesLabel -> {"Dx [mm]", "|kappa_|"}]

Export["KappaVsDx.eps", %]

dy[30] //. param

Plot[dy[dx] //. param, {dx, 0, 40}, 
      AxesLabel -> {"Dx [mm]", "Dy [mm]"}, 
      AxesOrigin -> {0, 0}, 
      PlotRange -> {{0, 40}, {-20, 30}}]

Export["DyDx.eps", %]

radius1[dxh_]  := Sqrt[(ay1 + dy[dxh])^2 + dxh^2]

radius2[dxh_]  := Sqrt[(ax2 -  dxh)^2 + dy[dxh]^2]

p10 = 
      Plot[radius1[dx] //. param, {dx, 0, 40}, 
        AxesLabel -> {"Dx [mm]", "Distance [mm]"}, 
        AxesOrigin -> {0, 0}, 
        PlotRange -> {{0, 40}, {0, 50}}]

p11  =  
      Plot[radius2[dx] //. param, {dx, 0, 40}, 
        AxesLabel -> {"Dx [mm]", "Distance [mm]"}, 
        AxesOrigin -> {0, 0}, PlotRange -> {{0, 40}, {0, 50}}, 
        PlotStyle -> Dashing[{0.02, 0.06}]]

Show[p10, p11]

Export["DyDxRadius.eps", %]

Table[{i , radius1[i] //. param, radius2[i] //. param}, {i, 5, 30}] // 
      TableForm

Table[{i , dy[i] //. param}, {i, 5, 30}] // TableForm

?Export

?Plot

BBLR1 beam sizes

sx1ave  =  (2.136 + 2.23)/2.

sy1ave  =  (2.302 + 2.207)/2.

sx2ave  =  (2.278 + 2.159)/2.

sy2ave  =  (2.374 + 2.068)/2.

dy2[dx_] := 
      dy //. FindRoot[((ay1 + dy)^2 - 
                      dx^2) ((ax2 - dx)^2 + 
                        dy^2)^2 + (dy^2 - (ax2 - 
                            dx)^2) (dx^2 + (ay1 + 
                              dy)^2)^2 == 0 //. 
            param, {dy, -50}]

Table[{i , dy2[i] //. param}, {i, 5, 30}] // TableForm

Plot[dy2[dx] //. param, {dx, 0, 40}, 
      AxesLabel -> {"Dx [mm]", "Dy [mm]"}, 
      AxesOrigin -> {5, 0}, 
      PlotRange -> {{5, 40}, {-50, 0}}]

Export["DyDx2.eps", %]

p5 = 
      Plot[Qx[dy2[dx], dx] //. param, {dx, 5, 30}, 
        PlotRange -> {{5, 30}, {- .0002, 0.0002}}]

p6  =  
      Plot[Qy[dy2[dx], dx] //. param, {dx, 5, 30}, 
        PlotRange -> {{5, 30}, {-0.0002, 0.0002}}]

Show[p1, p5]

dx1stexp  =  (ax2^2 - ay1^2)/(2 ax2)  //. param