شبیه‌سازی شبکه‌های بی سیم و کامپیوتری

شبیه ساز شبکه های کامپیوتری

https://netsimulate.net/simulation-and-analysis-of-computer-networks-in-matlab/

فیلم آموزشی در دست انجام

https://netsimulate.net/learning-simulate-wireless-networks-with-matlab-energy-approach/

مکان یابی اشیای درون ساختمان با WiFi و روش سه گانه در نرم افزار متلب

https://netsimulate.net/location-objects-inside-building-with-wifi-and-triple-method-in-matlab/

تشخیص گره های سیاه چاله در شبکه VANET با متلب

https://netsimulate.net/detection-black-hole-nodes-in-vanet-with-aodv-protocol-in-matlab/

و.

همواره مایلم با مشتریان خودم به شکل مستقیم ارتباط برقرار کنم، در همکاری های خودم با سایت ها معمولا این امکان برای من فراهم نمی‌شد به خاطر همین  ترجیح دادم که از ادامه همکاری با سایت های واسط خوداری کنم
برای من گرافیک زیبا، کدنویسی زیبا و ابشن دار(قابلیت تنظیم توسط کاربر) و خلاقیت اهمیت زیادی داره و در هر پروژه و شبیه سازی تمایل دارم که الگوریتمی بکر و نو خلق کنم.
یکی از ویژگی‌های اساسی کدنویسی بنده، نمایش کل سناریو شبیه‌سازی به شکل قدم به قدم و بصری هستش، در روزهای اینده فیلم‌های کوتاهی از برخی شبیه‌سازی اپلود می کنم تا بتونید بیشتر با این رویکرد بنده ارتباط برقرار کنید.

این دستورات خطای زمان اجرا را 100% از بین می‌برند اگر این دستورات بالای اسکریپت اصلی حودتون کپی و پیست کنین با هر بار اجرا حافظه و کنسول پاکسازی می شه و تمامی فرم‌ها و کادرها بسته می شند و مسیری که اسریپت اصلی از اونجا اجرا می شه به مراه تمام زیر پوشه‌ها در فهرست مسیرهای شناخته شده متلب اضافه می شه و  دیگه نیازی نیست مسیر به شکل دستی اضافه کنین.

ضمنا ای دستورات کاملا استاندارد نوشته شده و هم برای ورژن‌های قدیمی و هم ورژن جدید قابل استفاده است.

%% clear commands
delete(allchild(0));close all;
clear;
clc;

%% path commands
file=matlab.desktop.editor.getActive;
try 
fileDetail=dir(file.Filename);
folderName=fileDetail.folder;
catch
   index=strfind( file.Filename,'\');
   folderName=file.Filename(1:index(end)-1);
end
paths=genpath(folderName);
addpath(paths);
cd(folderName);

توضیحات تفضیلی در این رابطه را در این لینک بخوانید،یادآوری می‌شم سورس کدها به فروش نمی‌رسند(چون برای افراد یا سایت‌های مختلف نوشته شده.بحث کپی رایت) مگر با تغیرات و اعمال الگوریتم های جدید در ان‌ها.

classdef node
    
    properties
        x% x position
        y % y position
        z % z position
        energy % primary energy
        
    end
    
    
    methods
        
        function outObj=node(x,y,z,energy)
            outObj.x=x;
            outObj. y =y;
            outObj. z=z;
            outObj.energy =energy;
            
        end
        
        function out=getproperty(obj,in1)
            out.x=obj.x;
            out. y =obj.y;
            out.energy =obj.energy;
            
            
        end
    end
    
end

در کد بالا یک کلاس با عنوان نود ایجاد شده با ویژگی‌های موقعیت و سطح انرژِی، حالا می توان از این کلاس اشیایی ایجاد کرد و در یک یا چند متغیر یا در یک ارایه ذخیره کرد، به شکل زیر :
 

x=1;
y=3.4;
z=0;
energy=0.5;
net.nodes(NInd).nodObj=node(x2,y2,z2,energy);

در ارسال برودکست دو هدف وجود داره از طریق برودکست مقصد را پیدا کنیم، یا هدف ارسال همگانی باشه در شکل زیر مثالی از نوع اول را مشاهده می‌کنید که عموما در شبکه‌های بی‌سیم و خصوصا در الگوریتم مسیریابی AODV کاربرد داره.

 اگر سوالی دارید خوشحال می‌شم کمک کنم.

برای آموزش هر نوع شبکه‌‌ای که بتواند با تکیه بر ویژگی‌ها یک خروجی به ما بدهد، فقط همین یک جدول را نیاز داریم، به عبارتی برای آموزش هر نوع شبکه‌ای اعم از شبکه‌های عصبی، درخت تصمیم یا مثلا svm کافی است که یک جدولی مانند جدول بالا را داشته باشیم، سه ستون جدول بالا ویژگی‌ها و یک ستون به عنوان برچسب می‌باشد. مبنای این آموزش فقط به مشاهده توسط انسان یا نتایج آزمایش‌ها توسط انسان برمی‌گرده، البته این روش آموزش آموزش با نظارت است، اما نوع خاصی از آموزش را هم داریم که آموزش بدون نظارت است، در این حالت دیگه برچسب مشخص نمی‌شود و شبکه عصبی با توجه به روبط بین ویژگی‌ها و برخی قوانین معین شده، می‌تواند ورودی مثل سه ستون اول را دسته بندی کرده و نشان دهد که ویژگی‌های موجود در هر سطر مربوط به کدام گروه می‌‍باشد.

همواره مایلم با مشتریان خودم به شکل مستقیم ارتباط برقرار کنم، در همکاری های خودم با سایت ها معمولا این امکان برای من فراهم نمی شه.
برای من گرافیک زیبا، کدنویسی زیبا و ابشن دار(قابلیت تنظیم توسط کاربر) و خلاقیت اهمیت زیادی داره و در هر پروژه و شبیه سازی تمایل دارم که الگوریتمی بکر و نو خلق کنم.
یکی از ویژگی‌های اساسی کدنویسی بنده، نمایش کل سناریو شبیه‌سازی به شکل قدم به قدم و بصری هستش، در روزهای اینده فیلم‌های کوتاهی از برخی شبیه‌سازی اپلود می کنم تا بتونید بیشتر با این رویکرد بنده ارتباط برقرار کنید.

دلیل اصلی نوشتن این تابع، این بود که گاهی وقت‌ها لازم می شد که یک برداری از نظر تعداد تکرار عناصر،طول عناصر بدون تکرار و اندیس(موقعیت) تکرار عناصر مشخص بشه، که در متلب چنین تابعی نبود، اگر هم باشه من اطلاع ندارم، البته توابعی مثل intersect و ismember خروجی مشابهی دارند ولی این توابع  موقع مقایسه دو بردار می تونند چنین خروجی‌هایی را بدهند.

function [groupN,grValue,grReap,grindexList]=groupLable(trainLable )
if size(trainLable,1)>1
     trainLable=trainLable'; 
end
trainLable1=trainLable;
 
groupN=0;
while(~isempty(trainLable1))
    groupN=groupN+1;

    gr=(trainLable1==trainLable1(1));
    groupCluster(groupN).reap=sum(gr);
    groupCluster(groupN).grValue=trainLable1(1);


    trainLable1=trainLable1(~gr);
end
grValue=([groupCluster.grValue]);
grReap=([groupCluster.reap]);
grindexList=[];
for grL=grValue
   grindexList(end+1).grI=  (grL==trainLable).*(1:length(trainLable));
   grindexList(end).grI(grindexList(end).grI==0)= [];

end

end

نمونه خروجی و نحوه فراخوانی این تابع در تصویر زیر قابل مشاهده است:

یکی از توابعی که دو سال پیش توسعه دادم، تابع normalVect، این تابع  بسیاری کاربردی، یک ماتریس یک بعدی یا یک بردار را دریافت می کنه و مقادیر اون را به شکل پیش فرض بین 0 و 1 مپینگ( توزیع به زبان ساده‌تر) می‌کنه و اگر بازه را مشخص کنید بردار یا ماتریس را به همان بازه مپینگ می کنه. سه نمونه از کد فراخوانی تابع را در بالای اون تعبیه کردم ، که می تونید از داخل کنسول متلب یا داخل توابع یا اسکریپت ها او را فراخوانی کنید.
 

% give a vector or matrix for normalizing between tow number or 0 and 1
% vect=normalVect([1,-2,4,1,11,4,5])
% vect=normalVect([1,-2,4,1,11,4,5],-2,1)
% vect=normalVect([1,-2,4;11,4,5])

function vect=normalVect(varargin)
if nargin==0
    msgbox('give a 1 argoment atleast');
    vect=[];
elseif ~isnumeric(varargin{1})
    msgbox('abnormal argoment');
    vect=[];
else
   vect=varargin{1}; 
if nargin==1
    vect=(vect-min(vect(:)))./(max(vect(:))-min(vect(:)));
elseif nargin==3 ||  nargin>3
    if varargin{2}~=varargin{3}
        vect=(vect-min(vect(:)))./(max(vect(:))-min(vect(:)));
        vMin=min([varargin{2},varargin{3}]);
        vMax=max([varargin{2},varargin{3}]);
        vDif=vMax-vMin;
        vect=(vect*vDif)+vMin;
        if nargin>3
            msgbox('we used 3 argoment only');
        end
    else
        vect=(vect-min(vect(:)))./(max(vect(:))-min(vect(:)));
        msgbox('abnormal limits: normaled 0: 1');
    end
end

end
end

این دستورات خطای زمان اجرا را 100% از بین می‌برند اگر این دستورات بالای اسکریپت اصلی حودتون کپی و پیست کنین با هر بار اجرا حافظه و کنسول پاکسازی می شه و تمامی فرم‌ها و کادرها بسته می‌شند و مسیری که اسکریپت اصلی از اونجا اجرا می‌شه به همراه تمام زیر‌پوشه‌ها به فهرست مسیرهای شناخته شده متلب اضافه می شه و  دیگه نیازی نیست مسیرها را به شکل دستی اضافه کنین.

ضمنا این دستورات کاملا استاندارد نوشته شده، و هم برای ورژن‌های قدیمی و هم ورژن جدید قابل استفاده است.

%% clear commands
delete(allchild(0));close all;
clear;
clc;

%% path commands
file=matlab.desktop.editor.getActive;
try 
fileDetail=dir(file.Filename);
folderName=fileDetail.folder;
catch
   index=strfind( file.Filename,'\');
   folderName=file.Filename(1:index(end)-1);
end
paths=genpath(folderName);
addpath(paths);
cd(folderName);

classdef node
    
    properties
        x% x position
        y % y position
        z % z position
        energy % primary energy
        
    end
    
    
    methods
        
        function outObj=node(x,y,z,energy)
            outObj.x=x;
            outObj. y =y;
            outObj. z=z;
            outObj.energy =energy;
            
        end
        
        function out=getproperty(obj,in1)
            out.x=obj.x;
            out. y =obj.y;
            out.energy =obj.energy;
            
            
        end
    end
    
end

در کد بالا یک کلاس با عنوان نود ایجاد شده با ویژگی‌های موقعیت و سطح انرژِی، حالا می توان از این کلاس اشیایی ایجاد کرد و در یک یا چند متغیر یا در یک ارایه ذخیره کرد، به شکل زیر :
 

x=1;
y=3.4;
z=0;
energy=0.5;
net.nodes(NInd).nodObj=node(x2,y2,z2,energy);

NInd یک مقدار عددی صحیح مثبت، مثل 1، 2 و .

این تابع دارای چهار ورودی و چهار خروجی هست، ورودی اول  ایدی نودهای مقصد(مثلا سرخوشه) ورودی دوم ایدی نودهای مبدا(مثلا نودهای عادی که قرار از بین سرخوشه های یکی را انتخاب کنند، ورودی سوم فلق isplot( ایا مسیر بین مبدا و مقصدها رسم شود) ، ورودی چهارمisfigur ( ایا فرم جدیدی باز شود یا نه؟) خروجی اول unDestInd( مجموعه نودهای مقصد منتخب یونیک)، خروجی دوم gridDist( ماتریس فاصله بین نودهای مبدا و مقصد،  خروجی سوم selctedDest( بردار ایدی  مجموعه نودهای مقصد منتخب  غیر یونیک  ) خروجی چهارم selctedDict( کمترین فواصل انتخاب شده بین نودهای مبداء و مقصد).

%%caaling script
% find nearst dest point(source vector) from source point
%   sourcPos.x=randi([-100,100],1,100);
%   sourcPos.y=randi([-100,100],1,100);
 
% 
%   destPos.x=randi([-100,100],1,30);
%   destPos.y=randi([-100,100],1,30);
 
%   [unDestInd, gridDist,selctedDest,selctedDict]=getNearstDest2D(destPos, sourcPos,1,1);


%%


function  [unDestInd, gridDist,selctedDest,selctedDict]=getNearstDest2D(destPos, sourcPos,isplot,isFig)
gridDist=pdist2D(sourcPos,destPos  );%distance between all nodes with grid center points
[selctedDict,selctedDest]=min(gridDist,[],2);% here we used an simple algorithm to determine each node in network is belong to which cell in grid network
unDestInd=  unique(selctedDest);% if we unique the number of centers that are selected  by normal nodes, based on minimum distance, in reality we get the number coverage area
if nargin >=3
    if isplot==1
        if isFig
            figure;
        end
            plot( sourcPos.x, sourcPos.y,'ro') ;hold on
            plot( destPos.x ,destPos.y,'gd') ;hold on
            plot([sourcPos.x;destPos.x(selctedDest')] ,[sourcPos.y;destPos.y(selctedDest')],'k-');hold on;
    else
    end
end
end

 
% gridDist=pdist2D(sourcPos  );%distance between all nodes with destnation nodes
% gridDist=pdist2D(sourcPos,destPos  );%distance between all nodes within a vector
function gridDist=pdist2D(s, d   )

if nargin==1
    xTs=s.x' * ones(1,length(s.x) );
    yTs=s.y' * ones(1,length(s.y) );
 
    xTd=ones(1,length(s.x) )' * s.x  ;
    yTd=ones(1,length(s.y) )' * s.x  ;
 
    gridDist=sqrt((xTs - xTd ).^2 +(yTs - yTd ).^2 );
    
else
    
    
    xTs=s.x' * ones(1,length(d.x) );
    yTs=s.y' * ones(1,length(d.y) );

    xTd=ones(1,length(s.x) )' * d.x  ;
    yTd=ones(1,length(s.y) )' * d.y  ;

    gridDist=sqrt((xTs - xTd ).^2 +(yTs - yTd ).^2 );

end
end

نمونه خروجی این تابع

تو شبیه‌سازی شبکه‌های بی‌سیم ساخت اشیای سه بعدی و توزیع سه بعدی اهمیت زیادی داره، در فضای سه بعدی امکان شبیه‌سازی عملی و واقعی‌تر سیگنال و افت مسیر و چند مسیری و.  فراهم می‌شه،  در زیر کد ایجاد یک تایر به شکل سه بعدی در متلب آورده شده ، با کپی و پیست این کد در داخل یک اسکریپت متلب می‌توانید خروجی  زیر را مشادهده کنید.

%% clear commands
function xx
delete(allchild(0));close all;
clear;
clc;

%% path commands
file=matlab.desktop.editor.getActive;
try 
fileDetail=dir(file.Filename);
folderName=fileDetail.folder;
catch
   index=strfind( file.Filename,'\');
   folderName=file.Filename(1:index(end)-1);
end
paths=genpath(folderName);
addpath(paths);
cd(folderName);

 
side=10;
t=linspace(0,2*pi ,side+1);
t(6)=NaN;
side1=50;
t1=linspace(0, 2*pi ,side1+1);

  
 

 
% ph=0;
% x=r' *cos(t+ph);
% z=(0:side)' * ones(1,side+1);
% y=r' *sin(t+ph) ;


ph=0;
x= cos(t+ph) /5 ;
y=  0*(0:side)';  
z= sin(t+ph) /5  ;



% plot3(x,y,z);hold on

 
% path=(t1+ph).^2;path=path/max(path)

 path.x=[0,1,1,1,2,2,2,1,1,0];
 path.y=0:length(path.x)-1;
 
 path.x=  cos (t1+ph)*5 ;  
path.y= sin(t1+ph) *5  ;

 
    [path.x,path.y]=divider(path.x,path.y,2);
    [t1,t1]=divider(t1,t1,2);

%       plot(path.x ,path.y  , 'r');hold on
% xlabel( 'X'); ylabel( 'Y'); zlabel( 'Z')
 
     matrixObj.x= cos(t )  ;
    matrixObj.y=  0*(0:side)';
    matrixObj.z= sin(t  )   ;
    
for in=1:length(path.x)-1
    
    
    
    [ x, y, z]=rotateParameter(matrixObj,[0,0,1],rad2deg(t1(in)));
    
%      path.x=  cos (t + ph(in)) ;  
%      path.y=sin (t+ ph(in) )  ;
% 
     
%        plot3(x+path.x(in),y+path.y(in),z,'b-');hold on
       surf([x+path.x(in);x+path.x(in+1) ],[y+path.y(in);y+path.y(in+1) ],[z;z],'facecolor','k');hold on; alpha(0.2)
%       plot3(x ,y ,z,'b');hold on

    
   axis('equal') 
end

end

function [ newx, newy, newz]=rotateParameter(matrixObj,azel,alpha)

u = azel(:)/norm(azel);
alph = alpha*pi/180;
cosa = cos(alph);
sina = sin(alph);
vera = 1 - cosa;
x = u(1);
y = u(2);
z = u(3);
rot = [cosa+x^2*vera x*y*vera-z*sina x*z*vera+y*sina; .
    x*y*vera+z*sina cosa+y^2*vera y*z*vera-x*sina; .
    x*z*vera-y*sina y*z*vera+x*sina cosa+z^2*vera]';

x = matrixObj.x;
y = matrixObj.y;
z = matrixObj.z;

[m,n] = size(z);
if numel(x) < m*n
    [x,y] = meshgrid(x,y);
end

[m,n] = size(x);
newxyz = [x(:) , y(:) , z(:) ];
newxyz = newxyz*rot;
newx =   reshape(newxyz(:,1),m,n);
newy =  reshape(newxyz(:,2),m,n);
newz =  reshape(newxyz(:,3),m,n);

end

یکی از پرکاربردی ترین توابعی که تا حالا توسعه دادم تابع divider هستش، کار این تابع اضافه کردن تعداد معینی نقاط دو بعدی و سه بعدی روی خطوط دو بعدی و سه بعدی هستش، این تابع در دو مد کار می کند، حالت منصفانه که در آن بسته به طول هر خط روی خم، تعداد نقاط را روی آن اضافه می‌کنه، حالت غیر منصفانه که در آن روی تمام خطوط یک خم بدون توجه به طول آن به یک اندازه نقاط جدید اضافه می‌کنه

نمونه اجرای زیر حالت غیر منصافنه و منصفانه فراخوانی این تابع را نشان می‌دهد:

کد این تابع برای حالت دو بعدی

%%%%%%%%%%%%%%%%%%%%%%%%%%-----------software details-----------%%%%%%%%%%%%%%%%%%%%%%

%devloped in Iran(ardebil)                                       %
%Copyright (C) 2017(1396) by moussa hasanzadeh as the head devloper %%
%All rights reserved.                                             % 
%programed in mathlab 2016                                 
%gmail:mo30no@gmail.com
%phone 09147082079


%%%%%%%%%%%%%%%%%----------published by http://wcnet.blog.ir/------------%%%%%%%%%%%%%%


 



%version 1.1
%isFair=1 adding points is based on distance isFair=0 adding of points based on edges
function   [X,Y]=divider(X,Y,n,isFair)



try
    if isFair==1
    else
        isFair=0 ;
    end
catch
    
isFair=0;
end
switch isFair
    case 0
T1=0;
T2=0;

for i=1:length(X)-1
    temp1=0;
    temp1=linspace(X(i),X(i+1),n+2)  ;
    
    T1=[T1,temp1(1:end-1)];
end
T1=[T1,temp1(end)];
X=T1(2:end);


for i=1:length(Y)-1
    temp2=0;
    temp2=linspace(Y(i),Y(i+1),n+2)  ;
    
    T2=[T2,temp2(1:end-1)];
end
T2=[T2,temp2(end)];
Y=T2(2:end);
    case 1

dists= sqrt( (X(2:end)-X(1:end-1)).^2 + (Y(2:end)-Y(1:end-1)).^2) ;


n1=round((dists./min(dists))*n);
T1=0;
T2=0;

for i=1:length(X)-1
  
    temp1=0;
    temp1=linspace(X(i),X(i+1),n1(i)+2)  ;
    
    T1=[T1,temp1(1:end-1)];
end
T1=[T1,temp1(end)];
X=T1(2:end);


for i=1:length(Y)-1
    temp2=0;
    temp2=linspace(Y(i),Y(i+1),n1(i)+2)  ;
    
    T2=[T2,temp2(1:end-1)];
end
T2=[T2,temp2(end)];
Y=T2(2:end);
end


end


%%%%%%%%%%%%%%%%%%%%%%%%%%-----------software details-----------%%%%%%%%%%%%%%%%%%%%%%

%devloped in Iran(ardebil)                                       %
%Copyright (C) 2017(1396) by moussa hasanzadeh as the head devloper %%
%All rights reserved.                                             % 
%programed in mathlab 2016                                 
%gmail:mo30no@gmail.com
%phone 09147082079

% here  u can make net objects, setting them and connect them togather, send message and capture them ,. the last version of this library will release at the end of 1397(2018)
% the last version will be a combnation of real and virtual environment for Analysiing of network as a good ,.
 
%%%%%%%%%%%%%%%%%----------published by http://wcnet.blog.ir/------------%%%%%%%%%%%%%%

 

اسکریپت فراخوانی این تابع، برای تنظیم مد منصفانه و غیر منصفانه به جای isFair  به ترتییب مقدار 1 و 0 قرار دهید. مقدار n هم تعداد را مشخص می کنه

  net=[];
  net.width=100;
  net.length=100;
  net.side=7;
  net.nodNum=200;
  ph=pi/net.side;
  
  t=linspace(0,2,net.side+1);
  
  x=sin(pi*t+ph); x=(((x-min(x))/(max(x)-min(x)))-0.5)* net.length;
  y=-cos(pi*t+ph);y=(((y-min(y))/(max(y)-min(y)))-0.5)*net.width;

  
figure;plot(x,y,'-*');
[X,Y]=divider(x,y,5,1);
figure;plot(X,Y,'-*');






  net=[];
  net.width=100;
  net.length=100;
  net.side=7;
  net.nodNum=200;
  ph=pi/net.side;
  
  t=linspace(0,2,net.side+1);
  
  x=[1,3,1.2,6,4]
    y=[4,0,12.2,3,1]

  
figure;plot(x,y,'-*');
[X,Y]=divider(x,y,5,1);
figure;plot(X,Y,'-*');
axis('equal')

در کد زیر شما با اعمال تنظیمات می‌تونید تعداد نود(nodNum)، چند ضعلی بودن محیط شبکه (side)، رنگ بندی('r*') ، طول عرض محیط شبکه یا ناحیه شبکه(width و)  جهت گیری محیط شبکه(ph) را تغییر دهید، با این وجود اگر net.side را به 3، 5 و. تغییر بدهید، تعدادی از نودها در خارج از شبکه قرار می‌گیرند، برای حل این مشکل از تابع دیوایدر و یک الگوریتم ابتکاری می تونید استفاده کنید که در پست‌های بعدی در مورد ان شرح خواهم داد.

%% clear commands
delete(allchild(0));close all;
clear;
clc;

%% path commands
file=matlab.desktop.editor.getActive;
try 
fileDetail=dir(file.Filename);
folderName=fileDetail.folder;
catch
   index=strfind( file.Filename,'\');
   folderName=file.Filename(1:index(end)-1);
end
paths=genpath(folderName);
addpath(paths);
cd(folderName);


%% net building and initialization
  net=[];
  net.width=100;
  net.length=100;
  net.side=4;
  net.nodNum=200;
  ph=pi/net.side;
  
  t=linspace(0,2,net.side+1);
  
  x=sin(pi*t+ph); x=(((x-min(x))/(max(x)-min(x)))-0.5)* net.length;
  y=-cos(pi*t+ph);y=(((y-min(y))/(max(y)-min(y)))-0.5)*net.width;

  plot(x,y,'-');hold on;
  
  for NInd=1:net.nodNum
      x=(rand-0.5)*net.length;
      y=(rand-0.5)*net.width;
      energy=0.5;
      plot( x,y,'r*');hold on;
      
  end

این تابع کاربرد وسیعی در شبکه‌های بی‌سیم دارد،  و برای تشخیص وجود نقاط درون یا روی یک چند ضلعی مورد استفاده قرار می‌گیرد، در زیر یک تکه کد برای آشنایی بیشتر با این تابع به همراه توضیحات خط به خط آورده شده است.

 

side=18;     % side number of polygon
L = linspace(0,2.*pi,side+1); %Linearly spaced vector
xv = cos(L)'; % x aspect of polygon
yv = sin(L)';% y aspect of polygon
xq = randn(250,1);%x aspect of point vector
yq = randn(250,1);%y aspect of point vector
[in,on] = inpolygon(xq,yq,xv,yv);  % used for detection of Points located inside or on edge of polygonal region

figure

plot(xv,yv) % polygon
axis equal

hold on
plot(xq(in),yq(in),'r+') % points inside
plot(xq(~in),yq(~in),'bo') % points outside
hold off

خروجی اجرای این اسکریپت:

تو این فیلم یک خانه هوشمند شبیه‌سازی شده که در آن سه نفر در حال رفت و امد در محوطه و داخل خانه هستند، درها از طریق سنسور‌های حسگر حضور افراد باز و بسته می‌شند، هنگام ورود افراد به اتاق سنسورها حضور ان ها را حس کرده و به واحدی مرکزی اطلاع می‌دند  و واحد مرکزی بسته به شرایط  دستور روشن کردن چراغ اتاق را صادر می‌کند، در عین حال به شکل تصادفی از یک تا سه جا نشتی گاز رخ می‌ده یا ممکن اصلا نشتی رخ نده، ولی اگر نشتی داشته باشیم، سنسورهای حسگر آن را به واحد مرکزی گزارش می‌کنند، واحد مرکزی هم ابتدا گاز  و برق  را قطع کرده و سپس پنجره‌ها را روشن می‌کنند، در این شرایط حتی اگر افراد به داخل اتاق وارد شوند چراغ ها  روشن نمی شوند در زیر فیلم کوتاهی از ان شبیه‌سازی سه بعدی را می‌توانید مشاهده کنید.

نوشتن پشته پروتکلی، قانون خاصی  داره، دانستن این قانون ساده، خودش نصف مسیر!

برای این کار، باید ابتدا یک سرایند ایجاد کنید، سرایند در متلب اسم یک پوشه است که با علامت + شروع می‌شه، بعد در داخل این سرایند باید یک پوشه هم‌تام کلاس ایجاد کنید، سپس چند اسکریپت باز کرده و کد کلاس مربوط به هر لایه را در داخل اون میارید.

در پشته پروتکلی که خودم نوشتم،حدود 20 پروتکل اصلی اینترنت را در خودش داره، که از جمله می‌توان به پروتکل http, ip, icmp, tcp, udp ,. اشاره کرد، از یک سرایند تحت عنوان tcpIp استفاده کردم و چهار لایه کاربردی ، اینترنت، انتقال و لینک را در داخل اون پیاده‌سازی کردم.

شمای کله این پشته پروتکلی به شکل زیر:

 حالا در داخل هر یک از کلاس‌ها پروتکل‌های هر لایه در قالب توابع نوشته می‌شوند که برای نمونه می توان پروتکل dhcp را در دخل کلاس application به شکل زیر مشاهده کرد:
 

اما هر لایه یک تابع ورودی داره که معمولا به نام اون لایه inport اضافه می‌شه ، وظیفه این تابع بررسی اولیه پکت دریافتی از لحاظ نوع، محتوا و. هستش، سپس بر مبنای نتیجه حاصله، این تابع ، پکت ره به پروتکل‌های همان لایه یا به لایه بالاتر ارسال می‌کنه، زمانی که به لایه هدف رسید، که می تونه لایه لینک، اینترنت، انتقال یا هر لایه دیگه ای باشه، اون لایه بعد از ایجاد پاسخ اون را از همان مسیر به پایین ترین لایه ارسال می کنه تا در نهایت در اینترنت رها بشه و
ولی وقتی قرار یک کامپیوتر یا گوشی پکتی را ارسال کنه هر لایه از طریق تابعی هم نام لایه + outport استفاده می‌کنه تا در نهایت به پایین ترین لایه برسه در این حین از طریق لایه لینک در اینترت رها می شه

در ادامه اجرای پروتکل tcp, http , icmp, arp را با استفاده از شبیه‌سازه netGate_nbo تحت این پشته پروتکل، در فیلم کوتاه زیر مشاهده کنین، حجم این فیلم 9 مگابایت

یادآوری می‌شه که تمام شبیه‌سازی‌های مرتبط با شبکه کامپیوتر و امنیت در محیط این شبیه ساز که باز زبان متلب آن را توسعه دادم انجام می گیره.
برای سفارش پروژه های امنیت در شبکه‌های کامپیوتری می تونید با این شماره تماس بگیرید یا در واتس اپ و تلگرام پیام بدید 09147082079 

این درخت سه بعدی برای پروژه معماری یک دانشجوی ایرانی مقیم خارج  برای  یک سیستم ابرسانی پیاده سازی شده.
پیاده‌سازی با استفاده از یک تابع بازگشتی  و توابع سینوس و کسینوس انجام گرفته

برای اعمال بعد سوم از روش رشد مقطع استفاده شده، یعنی بعد از ایجا مسیر و مقطع به عنوان مثال دایره، مربع و مثلث،  این مقاطع در مسیر تعریف شده حرکت کرده و در نقاط دو بعدی یا سه بعدی کپی می‌شود.

با این وجود به خاطر بحث کپی رایت از انتشار کد در حال حاضر معذوریم.

یک راز مهمی که انمیشن های امروزی با کیفیت بالایی تولیدو توسعه داده می شوند، استفاده از حسگرهایی هست که حرکات  انسان و حیوان   حتی اندامکهای چهره و  را به کامپیوتر منتقل می کنند.
در این حالت حسگرها به بدن حیوان و انسان نصب می شوند و موقعیت فضایی بخش های مختلف بدن به یک مدل رایانه‌ای منتقل می شوند. و به عبارتی باهم دیگر ست می شوند، حالا با تغییر موقعیت انسان یا حیوان مدل رایانه ای رفتار و حرکات آنها را تقلید می‌کنه.

این حسگرها می توانند حسگرهای بی سیمی یا سیمی باشند  که کمی کند عمل می کنند یا حسگرهایی باشند که با استفاده از دوربین های مادونقرمز عمل می کنند و

نمونه ابتدائی ورود اطلاعات حرکات بدن در فیلم زیر را مشاهده کنید، همچنین کد متلب به همراه دیتاست پیاده روی و دویدن انسان می‌توانید دانلود کنید.

دو عدد اول از شما درخواست می کنه، بعد از آماده سازی کلید های خصوصی و عمومی یک متن را از شما خواهد خواست یا اینکه شما می تونید با یک سری تغییراتی یک بردار عددی را مستقیما به الگوریتم تحویل بدید، چون متن در نهایت به کد اسکی اون تبدیل می شه، بعد الگوریتم، متن رمزگذاری شده و متن رمزگشایی  شده را به همراه کد اسکی متن رمزگذاری شده و رمزگشایی شده به شما برمی گردونه.

clc;% clear theconsol
clear;%clear the memory and variables
delete(allchild(0)); %close all open diologs

file=matlab.desktop.editor.getActive;
try
      fileDetail=dir(file.Filename);
      folderName=fileDetail.folder;
catch
      index=strfind( file.Filename,'\');
      folderName=file.Filename(1:index(end)-1);
end
paths=genpath(folderName);
addpath(paths);
cd(folderName);

disp('RSA algorithm');
p=input('Enter the prime no. for p: ');
q=input('Enter the prime no. for q: ');
n=p*q;
fprintf('\nn=%d',n);
phi=(p-1)*(q-1);
fprintf('\nphi(%d) is %d',n,phi);

e=randsample( primes(n-1),1);

val1=0;
d=0;
while(val1~=1)
d=d+1;
val1=mod(d*e,phi);
end
fprintf('\nd=%d',d);
fprintf('\nPublic key is (%d,%d)',e,n);
fprintf('\nPrivate key is (%d,%d)',d,n);
m=input('\nEnter the message: ','s');
m1=m-0;
disp('ASCII equivalent of message ');
disp(m1);
over=length(m1);
 o=1;
 while(o<=over)
    m=m1(o);
    diff=0;
    if(m>n)
        diff=m-n+1;
    end
    m=m-diff;
      
    
    
    % beacause matlab dont support larg number modulus operation then
    % encription and deription is done in bit level
qm=dec2bin(e);
len=length(qm);
c=1;
xz=1;
while(xz<=len)  
    if(qm(xz)=='1')
      c=mod(mod((c^2),n)*m,n);
    elseif(qm(xz)=='0')
        c=(mod(c^2,n));
    end
    xz=xz+1;
end
c1(o)=c;
c1(o)=mod(m^e ,n);
% nm2(o)=mod(c1(o)^d ,n)+diff; 


% Public key is (e,n) => (7,33)
% Private key is (d,n) => (3,33)
% The encryption of m = 2 is c = 27 % 33 = 29
% The decryption of c = 29 id m =293 % 33 = 2




qm1=dec2bin(d);
len1=length(qm1);
nm=1;
xy=1;
while(xy<=len1)    
    if(qm1(xy)=='1')
       nm=mod(mod((nm^2),n)*c,n);
    elseif(qm1(xy)=='0')
        nm=(mod(nm^2,n));
    end
     xy=xy+1;    
end
nm=nm+diff;
nm1(o)=char(nm);
o=o+1;
 end
 
 
o=1;

fprintf('\nThe encrypted message in ASCII is \n');
 while(o<=over)
   fprintf('\t%d',c1(o)); 
   o=o+1;
end
o=1;

fprintf('\nThe encrypted message is \n');
disp(char(c1))


fprintf('\nThe decrypted message in ASCII is \n');
while(o<=over)
fprintf('\t%d',nm1(o));
o=o+1;
end
fprintf('\nThe decrypted message is: ');
disp(nm1);
fprintf('\n');

%%%%%%%%%%%%%%%%%%%%%%%%%%-----------simulation area-----------%%%%%%%%%%%%%%%%%%%%%%
% devloped in Iran(ardebil)                                 %
% Copyright (C) 2017(1396) by moussa hasanzadeh.           %
%All rights reserved.                                      %
%programed in mathlab 2016                                %
%gmail:mo30no@gmail.com
%phone: 09147082079


function vlia1
        clc
        clear;
        delete(allchild(0));
        sideN=4;
        t1=linspace(0,1,sideN+1);%vector between 0 and 1
        X=-sin(2*pi*t1+(pi/sideN));%node X coordinates with phaze shift
        Y=-cos(2*pi*t1+(pi/sideN));%node Y coordinates with phaze shift
        X=((((X-min(X))/(max(X)- min(X)))*2)-1)./5;%normalize  and schale X
        Y=((((Y-min(Y))/(max(Y)- min(Y)))*2)-1)./5;%normalize and schale Y
        Z= (1:length(Y))*0 ;

        for x=0:20:40
            for y=0:20:40
                for z=1:15
                    plot3(X+x,Y+y,Z+z,'k');hold on
                    plot3([X+x;X+x],[Y+y;Y+y],[Z+z*0;Z+z],'k');hold on
                end
            end
        end



        xVec=[0,20,40,40,20,0,0];
        yVec=[0,0,0,40,40,40,0];
        zVec1=yVec*0;
        zVec1([2,5])=zVec1([2,5])+3;
        [xVec,yVec]=dividerE(xVec,yVec,10);

        zVec=yVec*0;
        [~,zVec1]=dividerE(xVec,zVec1,10);



        hj1=fill3([xVec],[yVec],[zVec+5],'w');alpha(hj1 ,1)
        hj1=fill3([xVec],[yVec],[zVec+5+0.2],'k');alpha(hj1 ,0.2)

        hj3=fill3([xVec;xVec],[yVec;yVec],[zVec+5;zVec+5+0.2],'k');alpha(hj3 ,0.7)

        hj2=fill3([xVec],[yVec],[zVec+10],'w');alpha(hj2 ,1)
        hj2=fill3([xVec],[yVec],[zVec+10+0.2],'k');alpha(hj2 ,0.2)

        hj3=fill3([xVec;xVec],[yVec;yVec],[zVec+10;zVec+10+0.2],'k');alpha(hj3 ,0.7)

        hj3=fill3([xVec],[yVec],[zVec+15],'k');alpha(hj3 ,0.2)
        hj3=fill3([xVec;xVec],[yVec;yVec],[zVec+15;zVec+15+0.2],'k');alpha(hj3 ,0.7)

        xVecPentHous=[0,20,40,20,0];
        yVecPentHous=[0,0,0,0,0] ;
        zVecPentHous=[15,15,15,18,15];
        for mm=1:39
            % fill3(xVecPentHous,yVecPentHous+mm,zVecPentHous,'k--')
            km=fill3(xVecPentHous,yVecPentHous+mm,zVecPentHous,'g-');hold on
            alpha(km,0.2)
        end

        fill3([xVec],[yVec],[zVec+15],'r')



        zVec1=zVec1+15;

        fill3([xVec],[yVec],[zVec1],'r-','edgecolor','k')
        fill3([xVec],[yVec],[zVec1-0.2],'r-','edgecolor','k')

        hj1=fill3([xVec*2]-20,[yVec*2]-20,[zVec],'g');alpha(hj1 ,1)
        yardVect=[13:length(xVec),1:10]

        hj1=fill3([[xVec(yardVect)*2]-20;[xVec(yardVect)*2]-20],[[yVec(yardVect)*2]-20;[yVec(yardVect)*2]-20],[zVec(yardVect);zVec(yardVect)+2],'r');




        for difIndex=1:10
            zVec=zVec+0.2;
            hj1=plot3([xVec(yardVect)*2]-20,[yVec(yardVect)*2]-20,[zVec(yardVect)],'ks-','lineWidth',1);alpha(hj1 ,1);hold on

        end
        zVec(11:12)=zVec(11:12)+1;
        xVec(11)= xVec(10);
        xVec(12)= xVec(13);
       
        hj1=plot3([xVec*2]-20,[yVec*2]-20,[zVec],'k-','lineWidth',2);alpha(hj1 ,1);hold on
         zVec([10,13])=0;
        
         door=[10,11,12,13,10]
         matrixObj.x=[xVec(door)*2]-20;
         matrixObj.y=[yVec(door)*2]-20;
         matrixObj.z=[zVec(door)];
         [ newx, newy, newz]=rotateParameter(matrixObj,[0,0,1],90);
         
         hj1=fill3(newx-sqrt((xVec(10)-xVec(13)).^2)*1.3, newy-32, newz ,'w','lineWidth',2);alpha(hj1 ,1);hold on

        

        end


        function   [X,Y]=dividerE(X,Y,n)

        T1=0;
        T2=0;

        for i=1:length(X)-1
            temp1=0;
            temp1=linspace(X(i),X(i+1),n+2)  ;

            T1=[T1,temp1(1:end-1)];
        end
        T1=[T1,temp1(end)];
        X=T1(2:end);


        for i=1:length(Y)-1
            temp2=0;
            temp2=linspace(Y(i),Y(i+1),n+2)  ;

            T2=[T2,temp2(1:end-1)];
        end
        T2=[T2,temp2(end)];
        Y=T2(2:end);

        end






            function [ newx, newy, newz]=rotateParameter(matrixObj,azel,alpha)

            u = azel(:)/norm(azel);
            alph = alpha*pi/180;
            cosa = cos(alph);
            sina = sin(alph);
            vera = 1 - cosa;
            x = u(1);
            y = u(2);
            z = u(3);
            rot = [cosa+x^2*vera x*y*vera-z*sina x*z*vera+y*sina; .
                x*y*vera+z*sina cosa+y^2*vera y*z*vera-x*sina; .
                x*z*vera-y*sina y*z*vera+x*sina cosa+z^2*vera]';

            x = matrixObj.x;
            y = matrixObj.y;
            z = matrixObj.z;

            [m,n] = size(z);
            if numel(x) < m*n
                [x,y] = meshgrid(x,y);
            end

            [m,n] = size(x);
            newxyz = [x(:) , y(:) , z(:) ];
            newxyz = newxyz*rot;
            newx =   reshape(newxyz(:,1),m,n);
            newy =  reshape(newxyz(:,2),m,n);
            newz =  reshape(newxyz(:,3),m,n);

            end

این هم یک تصویر از نتیجه کار:

همانطور که می دونید متلب زوایه 360 درجه را به شکل دو بخش منفی و مثبت 180 درجه ای نشون می ده، این مسئله باعث بروز مشکل در پیاده سازی برخی الگوریتم‌هایی می شه که با زوایه کار می کنند، به عنوان مثال زمانی که شما تعدادی نقاط پخش شده در فضا دارید و بخواهید که کوتاهترین مسیر بین آنها را پیدا کنید، یک روش ساده مرتب سازی این نودها بر مبنای مختصات قطبی و خصوصا زوایه هستش، در چنین مواقعی نقاط قطبی دو بخشی پاسخگو نیست بر همین اساس من یک تابع ساده را برای حل این شکل نوشتم، که ورودی اون زوایه دو بخشی(0تا 180 و 0تا -180) و خروجی اون زوایه یک بخشی 0تا 360 درجه هستش.

کد متلب این تابع:

%convert 2 pole angle  to 1 pole degree or radian  
%angl=Angle2pol_To_1pol([5,95,-95,0],2)
%angl=Angle2pol_To_1pol( -95,2)
function angl=Angle2pol_To_1pol(angl,type)
if nargin==1
type=1;
end
switch type
    case 1% radian type
        angl=mod(angl,2*pi);
        miInd=angl<0;
        angl(miInd)=(2*pi)+angl(miInd);
        
    case 2
        angl=mod(angl,360);
        miInd=angl<0;
        angl(miInd)=(360)+angl(miInd);
end
end

نمونه استفاده از این تابع و خروجی حاصل در  تصویر زیر ابتدا مختصات دکارتی نکات درهم و برم به مختصات قطبی ان تبدیل شده و سپس با تبدیل زوایه دو بخشی به زاوایه یک بخشی 360 درجه و مرتب سازی آن بر اساس زاویه، شکل سمت چپ به شکل سمت راست تغییر پیدا کرده، این روش بسیار ساده می تونه برای پیدا کردن کوتاهترین مسیر هم مورد استفاده قرار بگیره.

همانطور که می‌دونید معادله عمومی منحنی بیضوی به شکل y² = x³ + ax + b می‌‍باشد که در ورژن بیت کوین مقدار ضریب a برابر صفر و مقدار ثابت b برابر 7 می‌باشد. نمای کلی منحنی بیضوی بیت کوین به شکل زیر می‌باشد:

از ویژگی‌های جالب این منحنی عبارتند از:
1-  اگر دو نقطه را در روی منحنی انتخاب کنید که خط واصل آن‌ها عمودی نباشد، خط واصل این دو نقطه، نقطه سومی را روی منحنی قطع خواهد کرد، که اگر این  نقطه را نسبت به محور xها متقارن کنیم، حاصل آن برابر با حاصل جمع دو نقطه اولیه خواهد بود.(point adding)

برای اضافه کردن دو نقطه به هم از روابط زیر استفاده می‌شود، که rx و ry حاصل اضافه کردن این دو نقطه می‌باشد.

c = (qy – py) / (qx – px)
rx = c2 – px – qx
ry = c (px – rx) – py

2- اگر یک نقطه ای را در روی این نمودار انتخاب کنیم، و چند بار ان نقطه را با خودش جمع کنیم، دوباره نقطه جدید روی این منحنی خواهد افتاد.(point doubling)

برای اضافه کردن یک نقطه به خودش از روابط زیر استفاده می‌شود، که rx و ry حاصل اضافه کردن این نقطه به خودش می‌باشد.

c = (3px2 + a) / 2py
rx = c2 – 2px
ry = c (px – rx) – py

در ادامه کد متلب پیاده‌سازی این عملیات را می‌توانید به شکل زیر مشاهده کنید:

clc;
clear;
close all;
% An elliptic curve is represented algebraically as an equation of the form
% y2 = x^3 + ax + b  ,For a = 0 and b = 7 (the version used by bitcoin), it looks like this y2 = x^3 + 7 

x=- nthroot(7,3):0.02: nthroot(7,3)*5; y = sqrt(((x.^3) +7));
 plot(x, real(y),'r- ', 'lineWidth',1);hold on;
 plot(x,-real(y),'r- ', 'lineWidth',1);hold on;



%% point adding
px=x(5) ;py=real(y(5)) ;qx=x(500) ; qy=real(y(500)) ;
c = (qy - py) /  (qx - px) ;
rx = (c^2)  - px - qx;
ry = c* (px-rx) -py;
 [X,Y]=polyxpoly( [px,qx],[py,qy],x,real(y)); 
plot([px,qx],[py,qy],'r*-')
plot( X,Y,'go-');hold on;
plot( rx,ry,'ko-');hold on;
text( rx,ry,'point adding');


%% point doubling
px=x(60) ;py=real(y(25)) ;
c = (3*(px^2) + 0) / (2*py);
rx = (c^2)  - (2*px);
ry = c*(px - rx) - py;

plot( px,py ,'cd-')
plot( rx,ry,'co-');
text( rx,ry,'point doubling');

نمونه خروجی:

این پست با پیاده‌سازی کامل بروزرسانی خواهد شد.

مطالعه بیشتر.

در صورت انتخاب موضوع و تایید استاد راهنما ، فورا به ما اطلاع بدین تا  وضعیت آن‌هارا بروزرسانی کنیم.

دسته 1: شبکه های بی سیم

1 مسیریابی در شبکه های ادهاک با استفاده از الگوریتم ،DSR+دایسترا 
2- مسیریابی در شبکه های ادهاک  با استفاده از الگوریتم AODV+دایسترا
3- مسیریابی در شبکه های ادهاک  با استفاده از الگوریتم ،HSR+دایسترا
4- مسیریابی در شبکه های ادهاک  با استفاده از الگوریتم ABR+دایسترا
5- مسیریابی در شبکه های ادهاک  با استفاده از الگوریتم FSR+دایسترا
6- مسیریابی در شبکه های ادهاک  با استفاده از الگوریتم WRP+دایسترا
7- مسیریابی در شبکه های ادهاک  با استفاده از الگوریتم CGSR +دایسترا
8- مسیریابی در شبکه های ادهاک  با استفاده از الگوریتم DSDV +دایسترا
9- مسیریابی در شبکه های ادهاک  با استفاده از الگوریتم GSR +دایسترا
10- مسیریابی در شبکه های ادهاک  با استفاده از الگوریتم TORA +دایسترا
11- مسیریابی در شبکه های ادهاک با استفاده از الگوریتم LMR +دایسترا

12-  مسیریابی در شبکه های حسگر بی سیم با استفاده از الگوریتم BEE (Balanced Energy-Efficient Clustering)
13- مسیریابی در شبکه های حسگر بی سیم با استفاده از الگوریتم BEEM (multi hop BEE) 
14-  مسیریابی در شبکه های حسگر بی سیم با استفاده از الگوریتم S-WEB (multi hop BEE) 
15-  مسیریابی در شبکه های بی سیم با استفاده از الگوریتم گرگ خاکستری
16-  مسیریابی در شبکه های بی سیم با استفاده از الگوریتم  CCS
17- مسیریابی در شبکه های بی سیم با استفاده از الگوریتم  EAR
18- مسیریابی در شبکه های بی سیم با استفاده از الگوریتم PGR
19- مسیریابی در شبکه های بی سیم با استفاده از الگوریتمPGR_EP(ترکیب دو الگوریتم بالا)

دسته 2: شبکه های کامپیوتری و امنیت

1- تشخیص نفوذ در شبکه های کامپیوتری با استفاده از  داده پردازی واکنشی
2- استخراج دانش از متن های رمز شده با استفاده از واقعیت کاوی داده
3-  کشف آلودگی شبکه های کامپیوتری براساس  نظریه بازی ها
4-  سیستم تشخیص نفوذ (IDS) چندگانه در شبکه های حسگر بی سیم بر اساس داده کاوی تکاملی

دسته 3: داده کاوی

1- بازسازی جملات یک متن در هم و برم با تکنیک داده پردازی
2- تشخیص اعتماد به نفس از روی چهره
3- تخمین ضریب هوشی  افراد از روی تحلیل  گنجینه لغات
4- تشخیص قدرت حافظه افراد از روی تحلیل جملات
5- تشخیص وضعیت روحی افراد از روی تحلیل سخن.
6- تشخیص وضعیت روحی افراد از روی چهره.
7- تشخیص احساس از روی چهره(تنفر، دوستی، همدردی و.)
8- تشخیص اضطراب ار وری چهره.
9- تشخیص نگرانی از روی تحلیل جملات.
 

دسته 4: اسنترنت اشیا

1-بهینه سازی مصرف انرژی با استفاده از اینترنت اشیا ر خانه های هوشمند

2- ایجاد پروتکل امن ارتباطی برای جلوگیری از هک دستگاه‌های الکترونیکی  در واحدهای هوشمند

ارائه یک سیستم تشخیص و منترل وضعیت بیمار در بیمارستان ها

3- استفاده از اینترنت اشیا برای تشخیص، و هره نگاری و تعقیب هوشمند افراد خرابکار در تظارها و راهپیمای ها

4- پیش بینی رفتار انسان براساس یادگیری عمیق

5- تشخیص میزان سلامتی انسان از روی پردازش چهره و حرکات بدن

6- تشخیص و دسته بندی نژادهای انسان با استفاده از پردازش چهره
 

دسته 5: این موضوعات به شکل دائمی  بروزرسانی خواهد شد

در صورت انتخاب موضوع و تایید استاد راهنما ، فورا به ما اطلاع بدین تا  وضعیت آن‌هارا بروزرسانی کنیم.
در صورتی که در حوزه کامپیوتر، نرم افزار و موضوعی نیاز دارین به شکل رایگان موضوعات برای شما پیشنهاد می کنیم، ارائه موضوعات از کارشناسی تا دکتر  قابل ارائه است.

شما می دونید که معادله یک خط در فضا به شکل y=ax+b هستش و  نمایش دیگر و کاربردی اون هم به شکل (y-y1=m(x-x1  هستش، حالا در اینجا ما فقط یک خط در دو بعد داریم، با این حساب ب داشتن شیب یا همان m  و یک نقطه (x1,y1) ما قادر خواهیم بود که با گرفتن مقدار مختصه x مقدار تابع را به ازای آن مختصه بدست بیاوریم. یا برعکس با داشتن مقدار y شما می‌توانید مقدار x  نقطه را به ازای آن بدست بیاوید . در زیز کد متلب این عمل را می آوریم.

%% clear commands
delete(allchild(0));% delete all figures and graphical objects
clear;% Remove items from workspace, freeing up system memory
clc;% clear the Console


m=5;
x1=5, y1=10;

syms x;
y=m*(x-x1)+y1;
x2=input('give a x coordinate:\n ');
x=x2;
y2=eval(y);
plot([x1,x2],[y1,y2],'o--g');hold on;
 
x3=input('give another one x coordinate:\n ');
x=x3;
y3=eval(y);
plot([x2,x3],[y2,y3],'o--m');

نمونه خروجی این کد:

طبق تصویر زیر برای تعیین اینکه نقطه در داخل مثلث قرار داره یا خارج از آن باید  زاویه سه مثلثی که رئوس مثلث با نقطه مورد نظر می‌سازند را محاسبه کنین، حالا براحتی  می‌تونین با جمع کردن زوایا، طبق شکل زیر موقعیت نقطه را مشخص کنین،  اگر α+β+γ=360 باشد، نقطه در داخل مثلث و اگر α+β+γ<360 در آنصورت نقطه در خارج از مثلث قرار دارد. از این ترفند ساده ریاضی برای تشخیص وجود نقاط در داخل چند ضلعی هم می‌توانید استفاده کنید.

کاربرد این مسئله ساده ریاضی بیشتر در چشم‌های الکترونیکی در ربات‌ها می‌باشد

کد متلب برای تصویر بالا بزودی منتشر خواهد شد.

The MATLAB code for the above image will be released soon.

طبق مثالی که اوردم ، سه نقطه به عنوان ورودی در نظر گرفته شده، که نقطه اول  C به عنوان مبداء و نقطه اتصال خط a و خط b ، نقطه B ( انتهای خط a) نقطه A، انتهای خط b. دقت کنین که زاویه‌ها با حروف بزرگ و خطوط روبروی هر زاوایه حرف کوچک آن زوایه می‌باشد.

شما با تغییر مقادیر نقاط می‌توانید زاویه‌های مرتبط با آن‌ها را مشاهده کنید

کد متلب این مثال:

%% clear commands
delete(allchild(0));% delete all figures and graphical objects
clear;% Remove items from workspace, freeing up system memory
clc;% clear the Console

%% path commands
file=matlab.desktop.editor.getActive;% get current script address
try % if dont occur any error this block will run
      fileDetail=dir(file.Filename);% get current script details
      folderName=fileDetail.folder;% get current script folder name
catch% else if an  error occurs then
      index=strfind( file.Filename,'\');% find index of back space(\) in path text
      folderName=file.Filename(1:index(end)-1);% select all path text from 1 to last Occurrence of \ as script folder name
end
paths=genpath(folderName);% make current path and all sub paths
addpath(paths);% add all paths in known paths of matlab
cd(folderName);% go to path that current script runned from that



 
A_x=5;
A_y=12;


B_x=A_x/2;
B_y=A_y*2;


C_x=0/3;
C_y=0;


c=sqrt( ((B_x-A_x)^2) +( (B_y-A_y)^2));
b=sqrt( ((A_x-C_x)^2) +( (A_y-C_y)^2));
a=sqrt( ((B_x-C_x)^2) +( (B_y-C_y)^2));
 tSlop=(a^2 +b^2 -c^2)/(2*a*b);
alpha_rad= (acos(tSlop));
alpha_deg=rad2deg(acos(tSlop));


figure;
plot([ C_x,A_x],[C_y,A_y],'r-');hold on
plot([ C_x,B_x],[C_y,B_y],'r-');hold on
plot([ A_x,B_x],[A_y,B_y],'b--');hold on

text(A_x,A_y,'A','Horizontalalignment', 'right');
text(B_x,B_y,'B','Horizontalalignment', 'left');
text(C_x,C_y,'C','Horizontalalignment', 'center');

text((A_x+C_x)/2 ,(A_y+C_y)/2,'b','Horizontalalignment', 'right');
text((B_x+C_x)/2 ,(B_y+C_y)/2,'a','Horizontalalignment', 'left');
text((A_x+B_x)/2 ,(A_y+B_y)/2,'c','Horizontalalignment', 'center');


text(C_x+2 ,C_y+2,[ '\alpha= ' num2str(alpha_deg)  '  deg '   '  =  ' num2str(alpha_rad)  '  rad '],'Horizontalalignment', 'center');

axis('equal')

خروجی نمونه(  output example )

برای اینکه مشخص شود که نقطه‌ای در زاویه دید یک نقطه مرجع دیگر با زاویه مشخص شده با دو خط قرار دارد یا نه، کافی که مختصات نقطه هدف و انتهای دو خط را مشخص  کرده و پس از انتقال مبداء به نقطه مرجع  و تبدیل آن‌ها به مختصات قطبی در 360 درجه، مقادیر قطبی انتهای دو خط با مقادیر نقطه هدف مقایسه شود.

 کد متلب کل عملیات فوق در زیر:(the related matlab code for this operation)

function fg=IsPointInsideAngle(basePoint,leftPoint,rightPoint,quPoint )

Px=basePoint.x;
Py=basePoint.y;


Ax=leftPoint.x;
Ay=leftPoint.y;

 Cx=rightPoint.x;
Cy=rightPoint.y;


 Dx=quPoint.x;
Dy=quPoint.y;

 
 

% transformed Ax and Cx
tAx= Ax-Px;
tAy= Ay-Py;

tCx= Cx-Px;
tCy= Cy-Py;

tDx= Dx-Px;
tDy= Dy-Py;
   
figure;plot([0, tAx],[0, tAy],'c');hold on;plot([0, tCx],[0, tCy],'c');hold on;plot(tDx,tDy,'rO');hold on
[tAtet,tArho]=cart2pol(tAx,tAy);tAtet=Angle2pol_To_1pol(tAtet,1);
[tCtet,tCrho]=cart2pol(tCx,tCy);  tCtet=Angle2pol_To_1pol(tCtet,1);
[tDtet,tDrho]=cart2pol(tDx,tDy); tDtet=Angle2pol_To_1pol(tDtet,1);

fg=isItBetween(tAtet,tDtet,tCtet);



end


function  fg=isItBetween(num1,quNum,num2)
fg=0;
[val,~]=sort([num1,num2],'ascend');
if  quNum>=val(1) && quNum<=val(2) 
      fg=1;
end
end

دقت کنین دو تابع بالا در یک اسکریپت هم نام تابع اول باید قرار گیرد وگرنه می تونید دو تابع را در دو اسکریپت مجزا هم نام توابع مذکور قرار دهید

برای فراخوانی این تابع می‌توانید از کد زیر استفاده کنید:( calling script of th above function)

basePoint.x=4;
basePoint.y=8;
          
leftPoint.x=16;
leftPoint.y=18;
          
rightPoint.x=20;
rightPoint.y=6;
          
quPoint.x=15;
quPoint.y=13;          
fg2=IsPointInsideAngle(basePoint,leftPoint,rightPoint,quPoint );

خروجی به شکل 1و 0 خواهد بود.

 the output will be as 1 , 0  or true and false

اغلب نواحی که برای شبیه سازی شبکه‌‌های بی‌سیم مورد استفاده قرار می‌گیرند، نواحی متقارن و عموما مربعی شکل می‌باشند ، حال آنکه به شکل طبیعی امکان توزیع متقارن نودها وجود نداره، در این پست سعی کردیم که کد توزیع نودها در فضاهای غیرمتقارن را منتشر کنیم، امیدواریم که براتون مفید واقع شده باشه.

در توزیع نودها از روش انتشار نقاط منحنی استفاده شده است.

کد متلب برای توضیحات بالا

%%%%%%%%%%%%%%%%%%%%%%%%%%-----------software details-----------%%%%%%%%%%%%%%%%%%%%%%
%devloped in Iran(ardebil)                                       
%Copyright (C) 2020(1398) by moussa hasanzadeh as the head devloper %%
%All rights reserved.                                            
%programed in mathlab 2016                                 
%gmail:mo30no@gmail.com
%phone 09147082079

%% clear commands
delete(allchild(0));% delete all figures and graphical objects
clear;% Remove items from workspace, freeing up system memory
clc;% clear the Console

%% path commands
file=matlab.desktop.editor.getActive;% get current script address
try % if dont occur any error this block will run
      fileDetail=dir(file.Filename);% get current script details
      folderName=fileDetail.folder;% get current script folder name
catch% else if an  error occurs then
      index=strfind( file.Filename,'\');% find index of back space(\) in path text
      folderName=file.Filename(1:index(end)-1);% select all path text from 1 to last Occurrence of \ as script folder name
end
paths=genpath(folderName);% make current path and all sub paths
addpath(paths);% add all paths in known paths of matlab
cd(folderName);% go to path that current script runned from that


%% net building  commands
net=[];% create an empty var as net object
net.para.width=1000;% width param of net
net.para.length=1000;% width param of net
net.para.height=5;% width param of net
net.para.side=6;% side of Polygon
net.para.nodNum=1000 ;% number of nodes of net work
net.para.simPer=1000;% number of simulation period
net.para.sInd=0;%period counter



 
ph=pi/net.para.side;% used for shifting of area to be Attractive!
t=linspace(0,1,net.para.side+1);% Generate linearly spaced vector( used for making polygon)
x1=sin(2*pi*t+ph); x1=(((x1-min(x1))/(max(x1)-min(x1)))-0.5)* net.para.length;% x coordinate of area(after normalising)
y1=-cos(2*pi*t+ph);y1=(((y1-min(y1))/(max(y1)-min(y1)))-0.5)*net.para.width;% y coordinate of area(after normalising)


  [x1,y1]=divider(x1,y1,10,1);% adding points between each pair of vertices

rs1=randsample(normalVect(rand(1,length(x1)),1,.5),length(x1));
figure;
x1=x1.*rs1;
y1=y1.*rs1;
[x1,y1]=divider(x1,y1,500);% adding points between each pair of vertices
x1(end+1)=x1(1);
y1(end+1)=y1(1);

plot(x1,y1,'r');hold on
rs=sqrt(randsample(normalVect(rand(1,length(x1)),.0,1),length(x1)));

plot(x1.*rs,y1.*rs,'.')






% give a vector or matrix for normalizing between tow number or 0 and 1
% vect=normalVect([1,-2,4,1,11,4,5])
% vect=normalVect([1,-2,4,1,11,4,5],-2,1)
% vect=normalVect([1,-2,4;11,4,5])

function vect=normalVect(varargin)
if nargin==0
    msgbox('give a 1 argoment atleast');
    vect=[];
elseif ~isnumeric(varargin{1})
    msgbox('abnormal argoment');
    vect=[];
else
   vect=varargin{1}; 
if nargin==1
    vect=(vect-min(vect(:)))./(max(vect(:))-min(vect(:)));
elseif nargin==3 ||  nargin>3
    if varargin{2}~=varargin{3}
        vect=(vect-min(vect(:)))./(max(vect(:))-min(vect(:)));
        vMin=min([varargin{2},varargin{3}]);
        vMax=max([varargin{2},varargin{3}]);
        vDif=vMax-vMin;
        vect=(vect*vDif)+vMin;
        if nargin>3
            msgbox('we used 3 argoment only');
        end
    else
        vect=(vect-min(vect(:)))./(max(vect(:))-min(vect(:)));
        msgbox('abnormal limits: normaled 0: 1');
    end
end

end
end





%version 1.1
%isFair=1 adding points is based on distance isFair=2 adding of points based on edges
function   [X,Y]=divider(X,Y,n,isFair)



try
    if isFair==1
    else
        isFair=0 ;
    end
catch
    
isFair=0;
end
switch isFair
    case 0
T1=0;
T2=0;

for i=1:length(X)-1
    temp1=0;
    temp1=linspace(X(i),X(i+1),n+2)  ;
    
    T1=[T1,temp1(1:end-1)];
end
T1=[T1,temp1(end)];
X=T1(2:end);


for i=1:length(Y)-1
    temp2=0;
    temp2=linspace(Y(i),Y(i+1),n+2)  ;
    
    T2=[T2,temp2(1:end-1)];
end
T2=[T2,temp2(end)];
Y=T2(2:end);
    case 1

dists= sqrt( (X(2:end)-X(1:end-1)).^2 + (Y(2:end)-Y(1:end-1)).^2) ;


n1=round((dists./min(dists))*n);
T1=0;
T2=0;

for i=1:length(X)-1
  
    temp1=0;
    temp1=linspace(X(i),X(i+1),n1(i)+2)  ;
    
    T1=[T1,temp1(1:end-1)];
end
T1=[T1,temp1(end)];
X=T1(2:end);


for i=1:length(Y)-1
    temp2=0;
    temp2=linspace(Y(i),Y(i+1),n1(i)+2)  ;
    
    T2=[T2,temp2(1:end-1)];
end
T2=[T2,temp2(end)];
Y=T2(2:end);
end


end

نمونه خروجی برای کد بالا(Output sample for the above code):

برای اینکه مشخص شود که نقطه‌ای در زاویه دید یک نقطه مرجع دیگر با زاویه مشخص شده با دو خط قرار دارد یا نه، کافی که مختصات نقطه هدف و انتهای دو خط را مشخص  کرده و پس از انتقال مبداء به نقطه مرجع  و تبدیل آن‌ها به مختصات قطبی در 360 درجه، مقادیر قطبی انتهای دو خط با مقادیر نقطه هدف مقایسه شود.

 کد متلب کل عملیات فوق در زیر:(the related matlab code for this operation)

function fg=IsPointInsideAngle(basePoint,leftPoint,rightPoint,quPoint )

Px=basePoint.x;
Py=basePoint.y;


Ax=leftPoint.x;
Ay=leftPoint.y;

 Cx=rightPoint.x;
Cy=rightPoint.y;


 Dx=quPoint.x;
Dy=quPoint.y;

 
 

% transformed Ax and Cx
tAx= Ax-Px;
tAy= Ay-Py;

tCx= Cx-Px;
tCy= Cy-Py;

tDx= Dx-Px;
tDy= Dy-Py;
   
figure;plot([0, tAx],[0, tAy],'c');hold on;plot([0, tCx],[0, tCy],'c');hold on;plot(tDx,tDy,'rO');hold on
[tAtet,tArho]=cart2pol(tAx,tAy);tAtet=Angle2pol_To_1pol(tAtet,1);
[tCtet,tCrho]=cart2pol(tCx,tCy);  tCtet=Angle2pol_To_1pol(tCtet,1);
[tDtet,tDrho]=cart2pol(tDx,tDy); tDtet=Angle2pol_To_1pol(tDtet,1);

fg=isItBetween(tAtet,tDtet,tCtet);



end


function  fg=isItBetween(num1,quNum,num2)
fg=0;
[val,~]=sort([num1,num2],'ascend');
if  quNum>=val(1) && quNum<=val(2) 
      fg=1;
end
end

دقت کنین دو تابع بالا در یک اسکریپت هم نام تابع اول باید قرار گیرد وگرنه می تونید دو تابع را در دو اسکریپت مجزا هم نام توابع مذکور قرار دهید

برای فراخوانی این تابع می‌توانید از کد زیر استفاده کنید:( calling script of th above function)

basePoint.x=4;
basePoint.y=8;
          
leftPoint.x=16;
leftPoint.y=18;
          
rightPoint.x=20;
rightPoint.y=6;
          
quPoint.x=15;
quPoint.y=13;          
fg2=IsPointInsideAngle(basePoint,leftPoint,rightPoint,quPoint );

خروجی به شکل 1و 0 خواهد بود.

 the output will be as 1 , 0  or true and false

 این کد نتیجه درستی نمی ده بزودی اصلاح و باز ارسال می شه.

this code hase incorrect result , it will correct and will send again.

برای اینکه مشخص شود که نقطه‌ای در زاویه دید یک نقطه مرجع دیگر با زاویه مشخص شده با دو خط قرار دارد یا نه، کافی که مختصات نقطه هدف و انتهای دو خط را مشخص  کرده و پس از انتقال مبداء به نقطه مرجع  و تبدیل آن‌ها به مختصات قطبی در 360 درجه، مقادیر قطبی انتهای دو خط با مقادیر نقطه هدف مقایسه شود.

 کد متلب کل عملیات فوق در زیر:(the related matlab code for this operation)

function fg=IsPointInsideAngle(basePoint,leftPoint,rightPoint,quPoint )

Px=basePoint.x;
Py=basePoint.y;


Ax=leftPoint.x;
Ay=leftPoint.y;

 Cx=rightPoint.x;
Cy=rightPoint.y;


 Dx=quPoint.x;
Dy=quPoint.y;
 

% transformed Ax and Cx
tAx= Ax-Px;
tAy= Ay-Py;

tCx= Cx-Px;
tCy= Cy-Py;

tDx= Dx-Px;
tDy= Dy-Py;


[tAtet,tArho]=cart2pol(tAx,tAy);tAtet=Angle2pol_To_1pol(tAtet,1);
[tCtet,tCrho]=cart2pol(tCx,tCy);  tCtet=Angle2pol_To_1pol(tCtet,1);
[tDtet,tDrho]=cart2pol(tDx,tDy); tDtet=Angle2pol_To_1pol(tDtet,1);

fg=isItBetweenAng(tAtet,tDtet,tCtet);

 
figure;plot([Px, Ax],[Py, Ay],'c');hold on;plot([Px,  Cx],[Py,  Cy],'c');hold on;plot( Dx, Dy,'rO'); 
 switch fg
 
       case 1
             title( 'it is in angle point'  )
       case 0
             title('it is  not in angle point')
 end

end



function  fg=isItBetweenAng(num1,quNum,num2)



fg=0;
[val,~]=sort([num1,num2],'ascend');

if (num1<=pi  && num2<=pi  ) || (num1 >=pi  && num2>=pi  )
      if  quNum>=val(1) && quNum<=val(2)
            fg=1;
      end
else
      
if val(2)-val(1)>pi
  if  (quNum<=val(1) && quNum>=0 ) ||   (quNum>=val(2) && quNum<=val(2) )
        fg=1;
  end
else
      if  (quNum>=val(1) && quNum<=pi ) ||   (quNum>=pi && quNum<=val(2) )
            fg=1;
      end
      
end
end
end

دقت کنین دو تابع بالا در یک اسکریپت هم نام تابع اول باید قرار گیرد وگرنه می تونید دو تابع را در دو اسکریپت مجزا هم نام توابع مذکور قرار دهید

برای فراخوانی این تابع می‌توانید از کد زیر استفاده کنید:( calling script of th above function)

basePoint.x=4;
basePoint.y=8;
          
leftPoint.x=16;
leftPoint.y=18;
          
rightPoint.x=20;
rightPoint.y=6;
          
quPoint.x=12;
quPoint.y=8;          
fg2=IsPointInsideAngle(basePoint,leftPoint,rightPoint,quPoint );

خروجی به شکل 1و 0 خواهد بود.

 the output will be as 1 , 0  or true and false

برداشت کلی این است که نوشتن یک تابع فیتنس عموما با الگوریتم‌های فراابتکاری معنی پیدا می‌کند، با این وجود یک تابع فیتنس چیزی نیست جز انتخاب جواب های بهینه از بین مجموعه‌هایی که به شکل تصادفی، حریصانه. انتخاب می‌شوند. تابع فیتنس می‌تواند از نوع مینیمم، ماکزیمم یا همگرا به سمت مقادیر یا بردارها و. باشد.

به عنوان یک مثال ساده برای تابع فیتنس برای ماکزیمم سازی  را  می‌توان به شکل زیر اورد.

 فقط حواستون باشه، برای نوشتن تابع فیتنس برای مینیمم سازی باید مقدار اولیه بردار پاسخ را برابر بی نهایت یا inf قرار بدید و برای ماکزیمم سازی این مقدار را برابر با 0 یا حداقل مقدار قرار بدید.

clc
clear;
close all;


[X,Y,Z] = peaks(100);
[val,ind]=min(Z(:));
v.val=val;
v.ind=ind;
cn=0;
hs=figure('name' ,  ['period: '  num2str(cn) ' of '  num2str((length(X(:))) )   ' curen  value:  '  num2str(v.val(end)),  ' real max value:  '  num2str( max(Z(:)))] );
 subplot(1,2,1);surf(X,Y,Z);hold on
 subplot(1,2,1);plot3(X(ind),Y(ind),Z(ind),'r*');hold on

while(true)
     ind= randi(length(X(:)),1);
     if  v.val(end)<Z(ind)
      v.val(end+1)=Z(ind);

          v.ind(end+1)= ind ;
 subplot(1,2,1);plot3(X(v.ind),Y(v.ind),Z(v.ind),'r*');hold on ; 
  subplot(1,2,2);plot3(X(v.ind),Y(v.ind),Z(v.ind),'r-*');hold on
  

     else
    disp('');
     end
      cn=cn+1;
 

     
 
     pause(0.002)

    if cn>=length(X(:)) || v.val(end)==max(Z(:))
           break;
    end
     hs.Name=['period: '  num2str(cn) ' of '  num2str((length(X(:))) )   ' curen  value:  '  num2str(v.val(end)),  ' real max value:  '  num2str( max(Z(:)))];
end

خروجی کد بالا به شکل زیر خواهد بود:

شما می توانید با تغییر معادله صفحه زیر  در سط 4،نمودار  سه بعدی  رشد یافته روی مسیر را روی آن صفحه مشاهده کنین.

f(p,q)=2*p  +2*q;

مثلا  می توانین  معادله صفحه در سط 4 را به فرم  2+ f(p,q)=2*p  +0.5*q;  تغییر دهین و نتیجه را مشاهده کنین.

حتی می تونین، به جای صفحه  معادله یک رویه را وارد کنین، مثلا سطر 4 را با معادله رویه f(p,q)=2*p  +.5*q+2*p*q; جایگزین کنین. و نتیجه شکل زیر را مشاهده کنین.

کد اجرایی:

clc;
clear;
close all;
syms p q X Y Z  yd(X) yl(X) f(p,q) ft(p,q) T xp(T) yp(T) zp(T)
f(p,q)=2*p  +2*q;
[p1,q1]=meshgrid( 0:.1:1 ,-1:.1:1 );
 surf(p1,q1,double(f(p1,q1)),FaceColor='g');alpha(.2);hold on;
 

t=linspace(0,1,100);
x=t;
y=sin(2*pi*x);yprim=2*pi*cos(2*pi*x);
plot3(x,y,f(x,y),'-r');hold on
n=length(t);
t1=linspace(0,1,200);
r=.08;
x1=r*sin(2*pi*t1);
y1=r*cos(2*pi*t1);
 
xM=[];
yM=[];
zM=[];

for s=1:length(t) -1

    %plot equLine for  slop and its 90 slop lines

    %% 
%     yprims =double((y(s+1)- y(s))/(x(s+1)- x(s)));
    yd(X)=yprim(s)*(X-x(s))+y(s);
    yl(X)=(-1./yprim(s))*(X-x(s))+y(s);
    tr=[x(s),x(s)+.00005 ];
    d1=diff( double(tr));
    d2=diff( double(yd(tr)));
    d3=diff( double(f(tr,yd(tr))));

    rd=[d1,d2,d3];
    xp(T)=x(s)+rd(1)*T;
    yp(T)=y(s)+rd(2)*T;
    zp(T)=double(f(x(s),y(s)))+rd(3)*T;
    tr1=double(solve((xp-x(s))^2+(yp-y(s))^2-r^2,T));
%     plot3(xp(tr1),yp(tr1) ,zp(tr1),'b-' );hold on
   


    %%

    d1=diff( double(tr));
    d2=diff( double(yl(tr)));
    d3=diff( double(f(tr,yl(tr))));

    r2=[d1,d2,d3];
%     r2=[1,double(yl(x(s)+1))-double(yl(x(s))),double(f(x(s)+1,double(yl(x(s)+1))))-double(f(x(s),double(yl(x(s)))))];
    xl(T)=x(s)+r2(1)*T;
    yl(T)=y(s)+r2(2)*T;
    zl(T)=double(f(x(s),y(s)))+r2(3)*T;
    tr2=(double(solve((xl-x(s))^2+(yl-y(s))^2-r^2,T)));
%     plot3(double(xl(tr2)),double(yl(tr2)) ,double(zl(tr2)),'c-' );hold on
    df=(double((yprim(s))))/abs(double((yprim(s))));
    
    if  df==-1
    tr2=tr2(end:df:1);
    end
    %%
%     plot3(x1+x(s),y1+y(s) ,double(f(x1+x(s),y1+y(s))),'b-' );hold on

    %%
    r3=double([subs(gradient(f ),[p,q],[x(s),y(s)])',-1]);
    xr=double(xl(tr2));
    yr=double(yl(tr2));
    zr=double(zl(tr2));

    xr0(T)=xr(1)+r3(1)*T;
    yr0(T)=yr(1)+r3(2)*T;
    zr0(T)=zr(1)+r3(3)*T; 
    
    xr1(T)=xr(2)+r3(1)*T;
    yr1(T)=yr(2)+r3(2)*T;
    zr1(T)=zr(2)+r3(3)*T; 
    tr3=[-.02,.02];
%     plot3(double([xr0(tr3),xr1([tr3(2),tr3(1)]),xr0(tr3(1))]),double([yr0(tr3),yr1([tr3(2),tr3(1)]),yr0(tr3(1))]),double([zr0(tr3),zr1([tr3(2),tr3(1)]),zr0(tr3(1))]));hold on
%     plot3(xr0(tr3),yr0(tr3),zr0(tr3),'*')
   
    xM=[xM;double([xr0(tr3),xr1([tr3(2),tr3(1)]),xr0(tr3(1))])];
    yM=[yM;double([yr0(tr3),yr1([tr3(2),tr3(1)]),yr0(tr3(1))])];
    zM=[zM;double([zr0(tr3),zr1([tr3(2),tr3(1)]),zr0(tr3(1))])];


  % peripendecular plane
  % peripendecular passpoints
    rpX=diff(double(xl(tr2))');
    rpY=diff(double(yl(tr2))') ;
    rpZ=diff(double(zl(tr2)'));
    B=[rpX,rpY,rpZ] ;
    
    C=double([subs(gradient(f ),[p,q],[x(s),y(s)])',-1]);
    C=cross(B,C);
    refPlane=C(1)*(X-x(s ))+C(2)*(Y-y(s))+C(3)*(Z-f(x(s ),y(s))) ;

    verPlaneCoff='XYZ';
    [~,ind]=find(C);
    ind=ind(end);
    in=find(~ismember([1,2,3],ind));

    ft(p, q)=subs(refPlane-C(ind)*sym(verPlaneCoff(ind)),[sym(verPlaneCoff(in(1))),sym(verPlaneCoff(in(2)))],[p,q]);
    ft=-(ft/C(ind));

      plot3(double([xr0(tr3),xr1([tr3(2),tr3(1)]),xr0(tr3(1))]),double([yr0(tr3),yr1([tr3(2),tr3(1)]),yr0(tr3(1))]),double([zr0(tr3),zr1([tr3(2),tr3(1)]),zr0(tr3(1))]),'r');hold on


% surf(double([xr0(tr3);xr1(tr3)]),double([yr0(tr3);yr1(tr3)]),double([zr0(tr3);zr1(tr3)]),FaceColor='g');hold on
% surf(double([xr0(tr3);xr1(tr3)]),double([yr0(tr3);yr1(tr3)]),double(ft(double([xr0(tr3);xr1(tr3)]),double([yr0(tr3);yr1(tr3)]))),FaceColor='b');hold on
plot3(double([xr0(tr3);xr1(tr3)])',double([yr0(tr3);yr1(tr3)])',double(ft(double([xr0(tr3);xr1(tr3)]),double([yr0(tr3);yr1(tr3)])))','b');hold on
plot3(double([xr0(tr3);xr1(tr3)]),double([yr0(tr3);yr1(tr3)]),double(ft(double([xr0(tr3);xr1(tr3)]),double([yr0(tr3);yr1(tr3)]))),'b');hold on



    view(3)
    xlabel('X');
    ylabel('Y');
    zlabel('Z');
%      axis equal
     pause(.0001)
end

axis equal
plot3(xM ,yM ,zM )
h=surf(xM ,yM ,zM ,'FaceColor','g');alpha(h,.1)