function $[$phases$,$ ref$\_$phase$\_$lags, phase$\_$lags, shil$\_$phase$\_$lags, binary$\_$states$]=$ SHIL$($model$,$ tspan, coupling, n$)$

T $=$ linspace$($tspan$(1),$ tspan$(2),1440)$;

dt $=$ mean$($diff$($T$))$;

$\%$ Reference with no coupling to compute period

sol $=$ model$(68,0,0,$ tspan, n$)$;

y $=$ sol.y;

t $=$ sol.x;

iy$1$ref $=$ interp$1($t$,$ y$(1,$:$),$ T$)$;

plot$($T$,$ iy$1$ref$)$;

$\%$ Compute period from last half of data

period $=$ compute$\_$period$($iy$1$ref$(721$:end$),$ dt$)$;

period

ts $=$ round$(1440-$ period $*6)$;

phase$\_$lags $=[]$;

ref$\_$phase$\_$lags $=[]$;

shil$\_$phase$\_$lags $=[]$;

binary$\_$states $=[]$; $\%$ To store binary states

$\%$ Input and reference signal

phases $=$ linspace$(0,$ pi$,8)$;

for phase $=$ phases

input$\_$signal $=$ prof$\_$pulse$($T$,$ period$/$n$,$ phase$*$n$)$;

ref$\_$signal $=$ prof$\_$pulse$($T$,$ period, phase$)$;

$\%$ Simulate system with input at SHIL frequency $($half natural period$)$

sol $=$ model$($period$,$ phase, coupling, tspan, n$)$;

y $=$ sol.y;

t $=$ sol.x;

iy$1=$ interp$1($t$,$ y$(1,$:$),$ T$)$;

figure$()$;

plot$($T$,$ iy$1)$;

$\%$ Compute phase lag of result to reference signal, using last half of data

phase$\_$lag $=$ compute$\_$phase$($iy$1($ts:end$),$ ref$\_$signal$($ts:end$),$ period, dt$)$;

phase$\_$lags$($end$+1)=$ phase$\_$lag;

$\%$ Compute phase lag of reference oscillator to ref signal, using last half of data

ref$\_$phase$\_$lag $=$ compute$\_$phase$($iy$1$ref$($ts:end$),$ ref$\_$signal$($ts:end$),$ period, dt$)$;

ref$\_$phase$\_$lags$($end$+1)=$ ref$\_$phase$\_$lag;

$\%$ Compute phase lag between SHIL and no SHIL

shil$\_$phase$\_$lag $=$ compute$\_$phase$($iy$1$ref$($ts:end$),$ iy$1($ts:end$),$ period, dt$)$;

shil$\_$phase$\_$lags$($end$+1)=$ shil$\_$phase$\_$lag;

$\%$ Determine binary state based on phase lag

threshold $=0$; $\%$ Adjust this threshold based on phase lag analysis

if phase$\_$lag $<$ threshold

binary$\_$states$($end$+1)=0$;

else

binary$\_$states$($end$+1)=1$;

end

figure$()$;

hold on;

plot$($T$($ts:end$),$ ref$\_$signal$($ts:end$)*$ max$($iy$1),\text{'}$r$--\text{'})$;

plot$($T$($ts:end$),$ input$\_$signal$($ts:end$)*$ max$($iy$1),\text{'}$g$--\text{'})$;

plot$($T$($ts:end$),$ iy$1($ts:end$),\text{'}$b$\text{'})$;

plot$($T$($ts:end$),$ iy$1$ref$($ts:end$),\text{'}$k$\text{'})$;

title$($gca$,$ sprintf$(\text{'}\%$f$\text{'},$ phase$\_$lag$))$;

end

ref$\_$

phase$\_$lags

phase$\_$lags

binary$\_$states $\%$ Display binary states

figure$()$;

plot$($ref$\_$phase$\_$lags, phase$\_$lags, $\text{'}$r$.\text{'},$ 'MarkerSize', $20)$;

hold on;

plot$($phases$*180/$pi$,$ phase$\_$lags, $\text{'}$g$.\text{'},$ 'MarkerSize', $20)$;

plot$([90,90],[0,180],\text{'}$k$--\text{'})$;

plot$([100,100],[0,180],\text{'}$k$--\text{'})$;

xlabel$(\text{'}$Initial phase $($degs$)\text{'})$;

ylabel$(\text{'}$Final phase $($degs$)\text{'})$;

end

can you give step by step explaination of this code