function [spos_n, svel_b, sCsn, simu]=la_imu(acc, gyro, gyro_der, pos_n, vel_b, Cbn, imu_la, imu_or, typ, gyro_vib, imu_vib)

nimu=size(imu_la,2);
spos_n=zeros(3,nimu);
svel_b=zeros(3,nimu);
sCsn=zeros(3,3,nimu);
simu=zeros(6,nimu);

%vibration information is not compulsory
if (isempty(gyro_vib))
    gyro_vib=zeros(3,nimu);
end

if (isempty(imu_vib))
    imu_vib=zeros(3,3,nimu);
    for in=1:nimu
        imu_vib(:,:,in)=eye(3);
    end
end

if (typ==0)
    [Rn, Re, g, sL, cL, WIE]=geo_param(pos_n);
end
for in=1:nimu
    %Attitude
    Crig=imu_or(:,:,in);
    Cvib=imu_vib(:,:,in);
    Csm=Crig*Cvib;
    Csn=Cbn*Csm;
    sCsn(:,:,in)=Csn;
    
    %position
    la_m=imu_la(:,in);
    if (typ==0)
        spos_n(:,in)=pos_n+[1/(Rn+pos_n(3)) 0 0;0 1/cL/(Re+pos_n(3)) 0; 0 0 -1]*Cbn*la_m;
    elseif (typ==1)
        spos_n(:,in)=pos_n+Cbn*la_m;
    end
    
    %Velocity
    if (typ==0)
        wie_m=Cbn'*[WIE*cL; 0; -WIE*sL];
    else
        wie_m=zeros(3,1);
    end
    svel_b(:,in)=Csm'*(vel_b+cross(gyro-wie_m,la_m));

    %Acceleration
    simu(1:3,in)=Csm'*(acc+cross(gyro,cross(gyro,la_m))+cross(gyro_der,la_m));

    %Rotation rate
    simu(4:6,in)=Csm'*(gyro+gyro_vib(:,in));
end