#############################
Script for making simulations
#############################

Simulations
-----------

::

   %Simulation of generic function (e.g. peaks) and S(Q,w) models

   %========= S(Q,w) model simulations ========================

   %simulate on sqw objects using S(Q,w) model for the physics

   %Use sqw_eval for S(Q,w) models
   sim_slice=sqw_eval(my_slice,@my_sqw_model,[1,2,3,4]);
   sim_slice_dnd=sqw_eval(d2d(my_slice),@my_sqw_model,[1,2,3,4]);

   %Here my_slice is an sqw object. We have a Matlab function on our path called my_sqw_model.m, and there are
   %some parameters associated with this model, [1,2,3,4].
   %
   %The form (1st line) of my_sqw_model.m should be:
   % function weight=my_sqw_model(h,k,l,e,pars)
   % the pars are the vector of input parameters we used above
   %
   % An example of an S(Q,w) model is given in the Horace "functions" directory, in the Matlab m-file
   % "demo_FM_spinwaves_2dSlice_sqw.m


   %Simulating a dnd and an sqw that are notionally the same gives different results, since for the sqw all of the contributing
   %detector pixels are simulated and then re-summed, whereas for the dnd just single points at the centre of bins are used
   %So if the intensity varies across a bin the dnd simulation will not capture this


   %========== Generic function simulation ==================

   %Use func_eval for generic functions
   peak_cut=func_eval(my_cut,@mgauss,[1,-1,0.1,2,0,0.1,1,1,0.1]);

   %Here the function mgauss is included in the Horace distribution (in the functions directory). It only works on 1-dimensional cuts, and simulates multiple
   %gaussians. The syntax is the same as for the sqw_eval function, so the vector argument is a list of input parameters for the function
   %In the case of mgauss, they are [amplitude1, centre1, width1, amplitude2, centre2, width2, ...., amplitudeM, centreM, widthM]
   %for M gaussians







Plotting dispersion
-------------------






::




   alatt=[2.87,2.87,2.87];
   angdeg=[90,90,90];

   lattice=[alatt,angdeg];
   rlp=[0,0,0; 1/2,0,0; 1/2,1/2,0; 1/2,1/2,1/2; 0,0,0; 0,0,1/2;];
   pars=[1,2,3,4,5];
   ecent=[0,0.1,10];
   fwhh=0.1;
   disp2sqw_plot(lattice,rlp,@my_dispersion_function,pars,ecent,fwhh);


   %The dispersion function has different (but similar) form to cross section functions
   %
   % The inputs are h,k,l and parameters; the outputs are the energy (i.e. dispersion) and spectral weight at that point

   %Here is an example dispersion function
   function [wdisp1,s_yy]=sr122_disp(qh,qk,ql,p)
   %
   % SrFe2As2 cross-section, from Tobyfit
   %

   %  Spin waves for FeAs, from Ewings et al., PRB 78
   %  Lattice parameters as for orthorhombic lattice i.e. a ~= b ~5.6Ang
   %
   % \tp(1)\tS_eff
   % \tp(2)\tSK_ab
   % \tp(3)\tSK_c
   %
   %    If ircoss=201:
   % \tp(4)\tSJ_1a
   % \tp(5)\tSJ_1b
   % \tp(6)\tSJ_2
   %
   %    If icross=202
   % \tp(4)\tS(2J2+J_1a)
   % \tp(5)\tS(2J2-J_1b)
   % \tp(6)\tSJ1a-SJ1b
   %
   % \tp(7)\tSJ_c
   % \tp(8)\tinverse lifetime gamma (= energy half-width)
   % \tp(19)\t0 if S(Q,w) as theory, =1 if multiply S(Q,w) by energy transfer
   % \tp(20)\t0 if twinned; 1 if twin #1 ; -1 if twin #2

   % If icross=207
   % As cross-section 204, but deal with J1a, J1b, J2 directly.

   s_eff = p(1);
   sj_1a = p(4);
   sj_1b = p(5);
   sj_2  = p(6);
   sj_c  = p(7);
   sjplus = sj_1a+(2.sj_2);
   sjminus = (2.sj_2)-sj_1b;
   sk_ab = 0.5.(sqrt((sjplus+sj_c).^2 + 10.5625) - (sjplus+sj_c));
   sk_c  = sk_ab;
   gam   = p(8);
   temp=4;

   alatt=[5.57,5.51,12.298];
   arlu=2pi./alatt;
   qsqr = (qh*arlu(1)).^2 + (qk*arlu(2)).^2 + (qlarlu(3)).^2;


   weight=zeros(size(qh));

   %First twin:
   a_q = 2.*( sj_1b.*(cos(pi.*qk)-1) + sj_1a + 2.*sj_2 + sj_c ) + (3.sk_ab+sk_c);
   d_q = 2.*( sj_1a.*cos(pi.*qh) + 2.*sj_2.*cos(pi.*qh).*cos(pi.*qk) + sj_c.*cos(pi.*ql) );
   c_anis = sk_ab-sk_c;

   wdisp1 = sqrt(abs(a_q.^2-(d_q+c_anis).^2));
   %wdisp2 = sqrt(abs(a_q.^2-(d_q-c_anis).^2));

   s_yy = s_eff.((a_q-d_q-c_anis)./wdisp1);