clear
set(0,'DefaultAxesFontSize',25)
set(0,'DefaultLinelinewidth',1)

% DESCRIPTION ::::::::::::::::::::
% A matlab routine to fit bolometric SN
% light curves with magnetar central engine model (with properties
% as in Kasen and Bildsten 2010) following the Arnett 1982
% formalism, with modification of kinetic energy treatment as described in
% Inserra, Smartt, Jerkstrand et al 2013.
% Author : A. Jerkstrand.
% Article reference : Inserra, Smartt, Jerkstrand et al ApJ 2013
% (denoted I13 below) 
% :::::::::::::::::::::::::::::::::

forceM='n';Mforce=0.1;
forcekappa='y';kappaforce=0.1;
forceE50='y';E50force=15;
forceslope='n';slopeforce=-1.4;

% CONTROLS ::::::::::::::::::::::::
obj = 'kilo'   % Select object to fit (4 char string).
alpha = 10;    % Power law exponent for envelope density profile
uselims = 0;   % 0: Ignore upper limits in fits. 1: Include upper limits in fits 
dotrap = 0;    % 0: Assume full X-ray/gamma trapping 1: Compute
               % trapping.
kappa_pp = 0.01; % Opacity to magnetar radiation (cm^2/g), only used if
                 % dotrap = 1
flow = 0.1;    % minimum multiplication factor from starting guess to probe
fhigh = 20;     % maximum multiplication factor from starting guess
               % to probe
N_Mej =30;     % Number of bins in Mej for coarse scan
N_E50 = 30;     % Number of bins in B14 for coarse scan
N_kappa = 30;     % Number of bins in Pms for coarse scan
Nslope = 30;    % Number of bins in time shift for coarse scan 

Nevals = 2000  % Parameter for fminsearch
tol = 1E-10    % Parameter for fminsearch

t_end = 20;   % Ending time of plotted light curve in days
save = 0;      % Save figures?
% -----------------------
% Usage : 
% Step 1 : Create a data input file along the following format:
%   First 2 lines : comments. Then
%   Column 1: Observer-frame epoch relative to peak (days)
%   Column 2: Logarithm of SN Rest-frame Bolometric luminosity (log10 L_bol (erg/s)) 
%   Column 3: Error on log10 L_bol
%   Any upper limits should be listed at the end, and controlled by
%   parameter "skip" below (=number of upper limits)
%
% Step 2: Copy and paste the section ("if (obj == 'test')") below
%   , changing the various variables to the new data set
%
% Step 3 : Set the fit variables above
%
% Step 4 : Run the code. Observed velocities and temperatures 
%   need to be added into Figs 2 and 3.
%
% Note that in the Arnett 1982 formalism, the light curve depends
% only on the power source and on the effective diffusion time
% scale tau_m = 1.05/sqrt(13.7*c)*kappa^(1/2)*Mej^(3/4)*Ek^(-1/4) s
% for constant density. The difference between the 56Ni/56Co case
% (as described in Arnett 1982) and the magnetar case is
% that energy input by the magnetar can be so large that the
% acceleration of the ejecta must be taken into account. We found
% by calibrating the solutions to the radiation hydrodynamic models
% in Kasen and Bildsten 2010 that using Ek = E0 + 1/2(Emag - Erad)
% gives excellent agreement : here E0 is the explosion energy, Emag
% is the total energy input by the magnetar, and Erad is the total
% radiated SN luminosity. In the fitting process here two steps are
% needed to get the final M_ej and velocities/temperatures correct:
% make first an initial fit with any value for E_51 (say 1) :
% obtain Emag and Erad from this run and plug into the expression
% for E_51 - rerun and the final M_ej value and
% temperature/velocity plots are obtained.

% Data and starting guesses for each SN:

%aa = fopen('/Users/anders/Dropbox/spectra_final/photometry_tables/PS17egl_bolom_dm33p11_red.txt')
%fgetl(aa); % Comments line 1
%fgetl(aa); % Comments line 1
%DATA = fscanf(aa,'%f',[8 Inf]); % Data, 3 columns (see above)

% ADD MORE DATA HERE AS FOR EXAMPLE ABOVE...

if (obj == 'kilo')

  flowM = flow
  fhighM = fhigh
  flowkappa= flow
  fhighkappa=fhigh
  flowslope = flow
  fhighslope = fhigh
  flowE50 = flow
  fhighE50 = fhigh

  %fid = fopen('/Users/anders/Dropbox/spectra_final/photometry_tables/PS17egl_bolom_dm33p11_red.txt'); % Input file
  fid = fopen('/Users/anders/Dropbox/spectra_final/photometry_tables/PS17egl_bolom_dm33p11_red_CI_v220170918-60547.txt')
  skip = 0;   % Number of upper limits (in end of input file) 
  redshift = 0.009; % redshift
  slopemin = -0.6
  slopemax = -2.5
  doptrap = 0
  if (dotrap == 0) % If assuming full trapping
      E50_kn = 1
      Mej_kn = 0.01
      kappa_kn = 1;
      slope_kn = -1.5
      savename = 'kilonova';

      if (forceM == 'y')
         Mej_kn = Mforce
	 flowM=1
	 fhighM=1
      end
      if (forcekappa == 'y')
         kappa_kn = kappaforce
	 flowkappa=1
	 fhighkappa=1
      end
      if (forceslope == 'y')
         slope_kn = slopeforce
	 flowslope=1
	 fhighslope=1
      end
      if (forceE50 == 'y')
         E50_kn = E50force
	 flowE50=1
	 fhighE50=1
      end 
      
  end 
  name = 'kilonova'
  save = 0;
  % For plotting:
  y0 = 40; % yrange-min
  y1 = 43; % yrange-max
  t1 = 42.8; t2= 42.5; t3=42.2; t4=42.9; % text positions
  xlv = 50;

end


% END OF USER INPUT ::::::::::::::::::::::::::::::::::::::::::::

forestgreen = [34 139 34]/255.;
gold = [0.75 0.75 0];

% READ IN DATA ------------------ :

fgetl(fid); % Comments line 1

fgetl(fid); % Comments line 2 

Z = fscanf(fid,'%f',[8 Inf]); % Data, 3 columns (see above)

fclose(fid);

tvec_days = Z(2,1:length(Z)-skip);  % Col 2 : Time in days observer frame
                                    % (arbitrary zero time),
                                    % Skip upper limits

ydata_log = Z(3,1:length(Z)-skip);  % Col 3 : Log10 L (erg/s)
                                    % bolometric
                                    % SN-frame lum. Skip
                                    % upper limits

sigma_log = Z(4,1:length(Z)-skip);  % Errors in log10 L_bol

N_data = length(tvec_days);  % Number of epochs

tvec_skip = [];

ydata_skip = [];

% ROUGH SPACE SCAN (exact best fit found by later optimization from this starting point) ------------------------

Mej_min = Mej_kn*flowM;  % Minimum ejecta mass
Mej_max = Mej_kn*fhighM;  % Maximum ejecta mass

E50_min = E50_kn*flowE50;  % Minimum B field in 10^14 G
E50_max = E50_kn*fhighE50; % Maximum B field in 10^14 G

kappa_min = kappa_kn*flowkappa;  % Minimun period in ms
kappa_max = kappa_kn*fhighkappa; % Maximum period in ms

slopemin = slope_kn*flowslope;
slopemax = slope_kn*fhighslope

% shiftmin and shiftmax are set above

Mej_vector = linspace(Mej_min, Mej_max, N_Mej);
E50_vector = linspace(E50_min, E50_max, N_E50);
kappa_vector = linspace(kappa_min, kappa_max, N_kappa);
slope_vector = linspace(slopemin, slopemax, Nslope);

bestchi2 = 1e50;  % Initialize best chi2 to a large value...

indb = find(tvec_days < 500); % TEMP

if (forceM == 'y')
   N_Mej =1
end
if (forcekappa == 'y')
   N_kappa =1
end
if (forceslope == 'y')
   Nslope =1
end
if (forceE50 == 'y')
   N_E50=1
 
end

for i = 1:N_Mej  % For each ejecta mass..

    for j = 1:N_E50  % For each energy...

        for k = 1:N_kappa   % For each spin period..

            for l = 1:Nslope   % For each time shift...

		%if (forceE50 == 'y')

		%  E50_vector(j) = 0.5*Mej_vector(i)*2e33*(0.2*3e10)^2/1e50
		    
		%end
  
                % Extract parameter values:
		x(1) = Mej_vector(i);
                x(2) = E50_vector(j);
                x(3) = kappa_vector(k);
                x(4) = slope_vector(l);

                vchar=sqrt(2*x(2)*1e50/(x(1)*2e33))/3e10;

	        if (vchar > 0.15 && vchar < 0.25) 

		  chi2(i,j,k,l) = chi2function_kilonova(x, kappa_pp, length(indb), tvec_days(indb), ydata_log(indb), sigma_log(indb), uselims, skip, tvec_skip, ydata_skip, redshift, dotrap);  % Compute the reduced chi2 value of this model
		 %pause
		else

		  chi2(i,j,k,l) = 1000000000000000.;

		end 
                                    
		if (chi2(i,j,k,l) < bestchi2)  % If new
                                                   % solution is
                                                   % the best one
                                                   % so far, store it...

                   bestchi2 = chi2(i,j,k,l);
                   bestMej = Mej_vector(i);
                   bestE50 = E50_vector(j);
                   bestkappa = kappa_vector(k);
                   bestslope = slope_vector(l);

                 end

            end

        end

    end

end

bestMej_coarse = bestMej;
bestE50_coarse = bestE50;
bestkappa_coarse = bestkappa;
bestslope_coarse = bestslope;
bestchi2_coarse = bestchi2;

disp('Best solution from coarse search: M_ej, B14, Pms, shift, chi2(reduced):')
bestMej
bestE50
bestkappa
bestslope
bestchi2
%pause

% Now perform an exact determination of the chi2 minimum: (cannot
% do globally from scratch, search often diverges);

x0 = [bestMej, bestE50, bestkappa, bestslope];  % Starting point is
                                              % best solution from
                                              % coarse search..

options = optimset('MaxFunEvals',Nevals);   % Increase maximum number of function evaluations, needed for some searches.
options = optimset('TolFun',tol);  % If upper limits are used,
                                   % this has to be small (the
                                   % default is 1e-4 and
                                   % optimization may terminate
                                   % prematurely for that value)



% Do minimization of chi2:

if (forceM == 'y')
   display('1')
   [x, fval] = fminsearchbnd(@(x) chi2function_kilonova(x, kappa_pp, length(indb), tvec_days(indb), ydata_log(indb), sigma_log(indb),uselims, skip, tvec_skip, ydata_skip, redshift, dotrap), x0, [Mforce 0.1 0.1 -3], [Mforce*1.01 10 100 -0.5], options) % Dont allow Pms below 1 ms

elseif (forcekappa == 'y' & forceslope == 'n' & forceE50 == 'y')
    display('x')
    %pause
    [x, fval] = fminsearchbnd(@(x) chi2function_kilonova(x, kappa_pp, length(indb), tvec_days(indb), ydata_log(indb), sigma_log(indb),uselims, skip, tvec_skip, ydata_skip, redshift, dotrap), x0, [1e-3 E50force kappaforce -3], [0.2 E50force*1.01 kappaforce*1.01 -0.5], options) % Dont allow Pms below 1 ms

elseif (forcekappa == 'y' & forceslope == 'y')
    display('aaa')
    %pause
    [x, fval] = fminsearchbnd(@(x) chi2function_kilonova(x, kappa_pp, length(indb), tvec_days(indb), ydata_log(indb), sigma_log(indb),uselims, skip, tvec_skip, ydata_skip, redshift, dotrap), x0, [1e-3 0.1 kappaforce slopeforce], [0.1 10 kappaforce*1.01 0.99*slopeforce], options)

elseif (forcekappa == 'y' & forceslope == 'n')
    display('bbb')
    %pause
    [x, fval] = fminsearchbnd(@(x) chi2function_kilonova(x, kappa_pp, length(indb), tvec_days(indb), ydata_log(indb), sigma_log(indb),uselims, skip, tvec_skip, ydata_skip, redshift, dotrap), x0, [1e-3 0.1 kappaforce -3], [0.1 10 kappaforce*1.01 -0.5], options)

elseif (forceslope == 'y' && forceE50 == 'n')
    display('apa')
    %pause
    [x, fval] = fminsearchbnd(@(x) chi2function_kilonova(x, kappa_pp, length(indb), tvec_days(indb), ydata_log(indb), sigma_log(indb),uselims, skip, tvec_skip, ydata_skip, redshift, dotrap), x0, [1e-3 0.1 0.1 slopeforce], [0.1 10 100 slopeforce*0.99], options) % Dont allow Pms below 1 ms

elseif (forceslope == 'y' && forceE50 == 'y')
  [x, fval] = fminsearchbnd(@(x) chi2function_kilonova(x, kappa_pp, length(indb), tvec_days(indb), ydata_log(indb), sigma_log(indb),uselims, skip, tvec_skip, ydata_skip, redshift, dotrap), x0, [1e-3 E50force 0.1 slopeforce], [0.1 E50force*1.01 100 slopeforce*0.99], options)

elseif (forceslope == 'n' && forceE50 == 'y')
  display('y')
  %pause
  [x, fval] = fminsearchbnd(@(x) chi2function_kilonova(x, kappa_pp, length(indb), tvec_days(indb), ydata_log(indb), sigma_log(indb),uselims, skip, tvec_skip, ydata_skip, redshift, dotrap), x0, [1e-3 E50force 0.1 -3], [0.2 E50force*1.01 100 -0.5], options)

else
   display('3')
   [x, fval] = fminsearchbnd(@(x) chi2function_kilonova(x, kappa_pp, length(indb), tvec_days(indb), ydata_log(indb), sigma_log(indb),uselims, skip, tvec_skip, ydata_skip, redshift, dotrap), x0, [1e-3 0.01 0.1 -3], [0.2 10 100 -0.3], options) % TEMP allowed limits
end
    			
if (fval < bestchi2) % Was a new best fit found?
 % Refined best fit values:
 bestMej = x(1);    
 bestE50 = x(2);
 bestkappa = x(3);
 bestslope = x(4);
 bestchi2 = fval;
 disp('Best solution exact search: M_ej, B14, Pms, shift, chi2(reduced):')
 bestMej
 bestE50
 bestkappa
 bestslope
 bestchi2
 Bestexact = x;

else

 disp('Refined search failed to improve solution.')

end
%pause
x = [bestMej, bestE50, bestkappa, bestslope];

%bestMej = 0.01; % TEMP
%bestE50 = 20; % TEMP
%bestkappa = 0.1 % TEMP
%bestslope=-1.2


% PLOTTING ----------

% Constants
c       = 2.998e10; % Speed of light (cm/s)   
e_ni0   = 3.90e10;  % Energy release in 56Ni decay (erg/s/g)  CHECKED. Note Arnett 1982 has too high value 4.8e10
e_co0   = 6.78e9;   % Energy release in 56Co decay (erg/s/g)  CHECKED 
m_sol   = 2e33;     % Solar mass (g) 
tau_ni  = 6.077*24*60*60/log(2); % Decay time (half life / ln2) for
                                 % Ni56 (s) 6.077 --> 8.8 d CHECKED
                                 % note ln in matlab is log

tvec = 10.^(linspace(-1,1.4))
				 
% Plot lightcurve ------------------------------

figure(1), clf, hold on, box on

  errorbar(tvec_days, ydata_log, sigma_log,'xr') % Convert time from observer-frame time to SN rest frame time

  % Plot upper limits, if there are any:
  if (skip > 0)

      %plot(tvec_skip/(1.+redshift), ydata_skip, 'rv','Markersize',10)
      plot(tvec_skip, ydata_skip, 'rv','Markersize',10)

  end

	%tt=linspace(1.4,8,300)
	%yy0 = log10(0.5e42*((tt/1.4).^(-0.6)))
	%yy = log10(0.5e42*((tt/1.4).^(-1.2)))
	%yy2 = log10(0.5e42*((tt/1.4).^(-2.4)))
	%plot(tt,yy0)
	%plot(tt,yy)
	%plot(tt,yy2)
	%tt=linspace(8,15,300)
	%yy0 = log10(0.4e41*((tt/8).^(-1.2)))
	%yy = log10(0.4e41*((tt/8).^(-2.5)))
	%yy2 = log10(0.4e41*((tt/8).^(-5)))
	%plot(tt,yy0)
	%plot(tt,yy)
	%plot(tt,yy2)
	
  %lvec=['0.001','0.01','0.1']
  stylevec = ['-','-','-',':',':',':']
  colorv = ['r','g','b','r','g','b']

  % Plot the best-fitting model:
  pls = '0'
  if (pls == '1')
   kappa=1
   x(1,1:4) = [0.001, 0.5, kappa, 0];
   x(2,1:4) = [0.01, 0.5, kappa, 0]; 
   x(3,1:4) = [0.1, 0.5, kappa, 0]; 
   x(4,1:4) = [0.001, 2, kappa, 0];
   x(5,1:4) = [0.01, 2, kappa, 0]; 
   x(6,1:4) = [0.1, 2, kappa, 0]; 
   Nmodels = 6
   for i = 1:Nmodels
     lcval(i,:) = lc_kilonova(x(i,:), M_ni, 1.0, kappa_pp, length(indb), tvec, dotrap)
   end
   for i = 1:Nmodels
     plot(tvec, log10(lcval(i,:)),'Linestyle',stylevec(i),'Color',colorv(i),'Linewidth',2)
     %set(h1,'linestyle',stylevec(i))
   end

   legend('Data','M=0.001','M=0.01','M=0.1')
   title('kappa=10')

  else

    bestMej,bestE50,bestkappa,bestslope, length(indb), tvec, dotrap
   lcval = lc_kilonova([bestMej,bestE50,bestkappa,bestslope],1.0,kappa_pp,length(indb),tvec,dotrap)

   plot(tvec, log10(lcval),'--')

  end

  tau_m = 1.05*((bestkappa/(13.7*c))^0.5) * (((((bestMej*m_sol)^3))/(bestE50*1e50))^0.25)
  %pause

  xlabel('Days since explosion');      

  ylabel('Log L (erg s^{-1})');

  %xmax = tend*1.1
  %xlim([-50 xmax])
  xlim([0.1 30])

  vchar = sqrt(2*bestE50*1e50/(bestMej*2e33))/3e10

  text(0.2,41.5,['M=',num2str(bestMej,2)],'Fontsize',20)
  text(0.2,41.0,['E50=',num2str(bestE50,2),'(v=',num2str(vchar,2),'c)'],'Fontsize',20)
  text(0.2,40.5,['kappa=',num2str(bestkappa,2)],'Fontsize',20)
  text(0.2,40.0,['slope=',num2str(bestslope,2)],'Fontsize',20)
  text(0.2,39.5,['chi2=',num2str(bestchi2,2)],'Fontsize',20)

  %text(6,41.5,'\beta=-0.6','Fontsize',20)
  %text(6,41.0,'\beta=-1.2','Fontsize',20)
  %text(4,40.3,'\beta=-2.4','Fontsize',20)

  %text(15,40.4,'\beta=-1.2','Fontsize',20)
  %text(15,40.0,'\beta=-2.5','Fontsize',20)
  %text(15,39.3,'\beta=-5','Fontsize',20)
    
  %ylim([y0 y1])
  ylim([39 43])

  set(gca,'Xminortick','on')

  set(gca,'Yminortick','on')

  set(gca,'Xscale','log')

  fs = 18

  xl = 145

  text(xl,t1,name,'Fontsize',fs+5,'Color','r')
  text(xl,t2,['M_{ej} = ' num2str(bestMej,2) ' M_{sun}'],'Fontsize',fs)
  %text(xl,t3,['B = ' num2str(bestB14,2) '*10^{14} G'],'Fontsize',fs)
  %text(xl,t4,['P = ' num2str(bestPms,2) ' ms'],'Fontsize',fs)

  if (save == 1)

      hgexport(gcf,savename, hgexport('factorystyle'),'Format','eps')

  end

% end of figure 1

% Compute photospheric velocities and temperatures:

Nmodels = 1
E_51(1) = bestE50*0.1
kappa = bestkappa

tthin_min = 10000000000
for j = 1:Nmodels

    Muse = bestMej
	
    V_ej = sqrt(10./3*E_51(j)*1e51/(Muse*2e33)); % Ejecta velocity, constant density  % CHECKED

    for i = 1:100 % For each time...

	      t = tvec(i)*24*3600;   % Time in seconds CHECKED

	      Rcore = V_ej*t;    % Outer ejecta radius in cm CHECKED

	      rhocore = bestMej*2e33/(4*pi/3.*Rcore^3);
	      %rhocore = x(j,1)*2e33/(4*pi/3.*Rcore^3);  % Core density in
                                               % g/cm^3 I13
                                               % Eq D12 CHECKED
 
	      taucore = bestkappa*Rcore*rhocore;   % Core optical depth in
                                     % optical CHECKED

	      tauatm = taucore/(alpha-1.);   % Atmosphere optical depth in
                                   % optical CHECKED

	      if (tauatm > 1)  % If atmosphere is still optically thick...

		  Rstar = Rcore*((alpha-1.)/taucore)^(1./(1.-alpha));  % Eq D19
                                                           % in I13
                                                           % CHECKED 

		  Tstar(j,i) = (lcval(i)/(4*pi*Rstar^2*5.67e-5))^0.25;  % Effective temperature CHECKED

		  Vstar(j,i) = Rstar/t;   % Photospheric velocity
				
	      elseif ((tauatm + taucore) > 0.5)   % Weird solutions when core is almost thin, unreliable solutions for such late phases

                  taumissing = 1 - tauatm;  % The optical depth that the core (between R* and Rcore) must provide

                  Rstar = Rcore - taumissing/(kappa*rhocore);  % I13 Eq D20 CHECKED

                  if (Rstar > 0.5*Rcore)   % Compute only if more than half the core is covered ny photosphere

		     Tstar(j,i) = (lcval(i)/(4*pi*Rstar^2*5.67e-5))^0.25;  % Blackbody formula CHECKED

		     Vstar(j,i) = Rstar/t; % CHECKED

                  else

	             Tstar(j,i) = 0;  % Stop when more than half core has been exposed

                     Vstar(j,i) = 0;

                  end

              else  % Weird solutions when core is almost thin, unreliable solutions for such late phases

		  if (t < tthin_min)

		     tthin_min = t/(24*3600)

		  end 

		  Tstar(j,i) = 0;

                  Vstar(j,i) = 0;

             end

      end

end

figure(1), hold on
plot([tthin_min,tthin_min],[39,43],'k--')
plot([tau_m/(24*3600), tau_m/(24*3600)],[39,43],'b--')

% Write to file

fid = fopen('output','w')

 for i = 1:length(tvec)

   fprintf(fid,'%12.8e %12.8e %12.8e %12.8e\n', tvec(i), lcval(i), Vstar(1,i), Tstar(1,i))

 end

fclose(fid)

% PLOT EFFECTIVE TEMPERATURES -------------------------

figure(2), clf, hold on

  for i = 1:Nmodels

    ind = find(Tstar(i,:) > 0);

    plot(tvec(ind), Tstar(i,ind))

  end

  errorbar(tvec_days,ydata_log,sigma_log,'*')

  set(gca,'Yminortick','on')

  title('T_{eff}')

  % Add plotting of observed temperatures here... 

  xlabel('Days since explosion')

  set(gca,'Xminortick','on')

  set(gca,'Xscale','log')

  ylim([0 2e4])
  xlim([0.1 30])

  if (save == 1)

    hgexport(gcf,strcat(savename,'_T'), hgexport('factorystyle'),'Format','eps')

  end

% PLOT PHOTOSPHERIC VELOCITIES----------------------

figure(3), clf, hold on, box on

  plot([0.1,20],[30,30],'--')
  text(10,35,'0.1c','Fontsize',20)

  for i =1:Nmodels
      
     ind = find(Vstar(i,:) > 0);

     plot(tvec(ind), Vstar(i,ind)/1e5/1e3)

  end

  %plot(tvec,3*ones(length(tvec)),'--')
  
  fs = 15
  
  %errorbar(timetemp, vfeii/1e4,errorfeii/1e4, 'ok','Markersize',fs,'Markeredgecolor','k','Markerfacecolor',forestgreen)
  %errorbar(timetemp, vmgii/1e4,errormgii/1e4, 'xk','Markersize',fs,'Markerfacecolor',forestgreen,'Markeredgecolor',forestgreen)
  %errorbar(timetemp, voi/1e4,erroroi/1e4, 'dk','Markersize',fs,'Markerfacecolor',forestgreen,'Markeredgecolor','k')
  %fs = 15;
  %text(xlv,1.7e4/1e4,name,'Fontsize',30,'Color',forestgreen)

  xlim([0.1 30])

  set(gca,'Xscale','log')
  
  xlabel('Rest frame days')

  ylabel('V ($10^3$ km s$^{-1})$','Interpreter','latex')

  set(gca,'Xminortick','on')

  set(gca,'Yminortick','on')

  ylim([0 100])

  if (save == 1)

       hgexport(gcf,strcat(savename,'_v'), hgexport('factorystyle'),'Format','eps')

  end

      
