Scilab – Uma abordagem prática aplicada a problemas reais da engenharia.

O objetivo deste trabalho é apresentar a aplicação da metodologia de aprendizagem baseada em problemas (PBL) para aulas da disciplina de Algoritmos e Técnicas de Programação em Scilab, para cursos de Engenharia, com o intuito de comprometer os alunos com a resolução de problemas reais de engenharia, através do uso da PBL de tal forma que os mesmos sintam-se inspirados a participar das aulas e desta forma, possam ter um aproveitamento melhor em contraste aos métodos tradicionais de ensino e aprendizagem baseada na criação de programas de proposito geral em linguagem de programação.

O livro está organizado em 4 capítulos, sendo o capítulo 1 uma introdução, explicando a origem e a história do Scilab e como ele se tornou essa excelente ferramenta de programação para solução de problemas numéricos.

No capítulo 2 são abordadas as fundamentações teóricas necessárias para o conhecimento da PBL e como ela pode ser útil na solução de problemas de engenharia. Ainda neste capítulo, são apresentados os métodos sobre elementos finitos, diferença finitas e métodos de iteração.

O capítulo 3 é reservado ao aprendizado do ambiente do Scilab, onde os principais comandos são explicados de forma clara. Vários pequenos exemplos são colocados para exemplificar os comandos explicados.

No capítulo 4, diversos problemas da área de engenharia são abordados, todos com suas devidas fundamentações matemáticas e proposta de soluções com elaboração de programas no Scilab, apresentação e discussão dos resultados.

Os programas descritos no livro, podem ser baixados através desse link.

Adicionalmente, os programas utilizados no livro, estão disponíveis logo abaixo, para que o leitor possa copiá-los e colá-los em seu programa do Scilab instalado no seu computador.

O livro está disponível para venda no site do clube dos autores e pode ser adquirido através deste link.


Página 77 – Efeito de transmissão térmica

clear;
clc;

t1=input("Qual é a temperatura inicial?"); //34
t2=input("Qual é a temperatura final?"); //60
A=[t1 0 0 0 0 0 t2]; // Declara o vetor A
plot2d(A);//Gráfico onde Y é temperatura e x o comprimento
aguarde=input("Iteração 0 – Digite ENTER"); //Pausa

// Realiza a primeira ITERAÇÃO
clc;
A(2)=(A(1)+A(3))/2; // Calcula o potencial de A(2)
A(3)=(A(2)+A(4))/2; // Calcula o potencial de A(3)
A(4)=(A(3)+A(5))/2; // Calcula o potencial de A(4)
A(5)=(A(4)+A(6))/2; // Calcula o potencial de A(5)
A(6)=(A(5)+A(7))/2; // Calcula o potencial de A(6)
plot2d(A);//Gráfico onde Y é temperatura e x o comprimento
aguarde=input("Iteração 1 – Digite ENTER"); //Pausa

// Realiza a segunda ITERAÇÃO
clc;
A(2)=(A(1)+A(3))/2; // Calcula o potencial de A(2)
A(3)=(A(2)+A(4))/2; // Calcula o potencial de A(3)
A(4)=(A(3)+A(5))/2; // Calcula o potencial de A(4)
A(5)=(A(4)+A(6))/2; // Calcula o potencial de A(5)
A(6)=(A(5)+A(7))/2; // Calcula o potencial de A(6)
plot2d(A);//Gráfico onde Y é temperatura e x o comprimento
aguarde=input("Iteração 2 – Digite ENTER"); //Pausa

// Realiza a terceira ITERAÇÃO
clc;
A(2)=(A(1)+A(3))/2; // Calcula o potencial de A(2)
A(3)=(A(2)+A(4))/2; // Calcula o potencial de A(3)
A(4)=(A(3)+A(5))/2; // Calcula o potencial de A(4)
A(5)=(A(4)+A(6))/2; // Calcula o potencial de A(5)
A(6)=(A(5)+A(7))/2; // Calcula o potencial de A(6)
plot2d(A);//Gráfico onde Y é temperatura e x o comprimento
aguarde=input("Iteração 3 – Digite ENTER"); //Pausa

// Realiza a quarta ITERAÇÃO
clc;
A(2)=(A(1)+A(3))/2; // Calcula o potencial de A(2)
A(3)=(A(2)+A(4))/2; // Calcula o potencial de A(3)
A(4)=(A(3)+A(5))/2; // Calcula o potencial de A(4)
A(5)=(A(4)+A(6))/2; // Calcula o potencial de A(5)
A(6)=(A(5)+A(7))/2; // Calcula o potencial de A(6)
plot2d(A);//Gráfico onde Y é temperatura e x o comprimento
aguarde=input("Iteração 4 – Digite ENTER"); //Pausa

Página 88 – Distribuição de potenciais ao longo de um condutor

clear; // Limpa a memória
clc; // Limpa o console
t1=input("Qual é a temperatura inicial?"); //34
t2=input("Qual é a temperatura final?"); //60
n=input("Qual é o número de elementos?"); //22
iteracao= input("Quantas épocas"); // Iterações

// Estrutura de repetição para preparar o vetor A
// Inicia em 1, incrementa 1 a cada repetição e
// termina em n.
for j=1:1:n // j variável contadora
	A(j)=0; // Zera cada elemento de A
end // Será repetido até aqui
A(1)=t1; // Atribui o valor da temperatura t1
A(n)=t2; // Atribui o valor da temperatura t2
plot2d(A);//Gráfico onde Y é temperatura e x o comprimento
aguarde=input("Iteração 0 – Digite ENTER"); //Pausa

// Realiza todas as iterações
clc;
for j=1:1:iteracao // vai da 1° a última iteração
	for i=2:1:(n-1) // vai do 2° ao elemento n-1
		A(i)=( A(i-1)+ A(i+1))/2; // Met. Iteração
	end // Fim do laço de repetição i
	plot2d(A); //Gráfico onde Y é temp. x o comprimento
	mprintf("Epoca=%d, ",j); //Informa a iteração corrente
	aguarde=input("Iteração 0 – Digite ENTER"); //Pausa
end // Fim do laço de repetição j
mprintf("Bom dia \n Fim da analise..."); // fim

Página 100 -Busca por solução sub-ótima

clear; // Limpa a memória
clc; // Limpa o console
t1=input("Qual é a temperatura inicial?"); //34
t2=input("Qual é a temperatura final?"); //70
n=input("Qual é o número de elementos?"); //50
iteracao= input("Quantas épocas"); // 10000
T= input("Qual é a precisão (0.1 a 0.001) // Precisão

// Estrutura de repetição para preparar o vetor A
// Inicia em 1, incrementa 1 a cada repetição e
// termina em n.
for j=1:1:n // j variável contadora
	A(j)=0; // Zera cada elemento de A
end // Será repetido até aqui

A(1)=t1; // Atribui o valor da temperatura t1
A(n)=t2; // Atribui o valor da temperatura t2
plot2d(A);//Gráfico onde Y é temperatura e x o comprimento
aguarde=input("Iteração 0 – Digite ENTER"); //Pausa
// Realiza todas as iterações
clc;
for j=1:1:iteracao // vai da 1° a última iteração
	Dif=0;
	for i=2:1:(n-1) // vai do 2° ao elemento n-1
		Aant(i)=A(i); // Grava cada ponto antes da mudança
		A(i)=(A(i-1)+A(i+1))/2; // Met. Iteração
		Dif=Dif+(A(i)-Aant(i)); // Soma das diferenças
	end // Fim do laço de repetição i
	Dif=Dif/n;
		if Dif<T then // Executa se a condição for atendida
		plot2d(A); //Gráfico onde Y é temp. x o comprimento
		mprintf("\nEpoca=%d, Diferença=%f",j,Dif);
		mprintf("\nEncontrada a solução sub-ótima! \n");
		break; //Interrompe o programa pois achou a solução
	end // fim do condicional
	plot2d(A); //Gráfico onde Y é temp. x o comprimento
	clc;
	mprintf("Epoca=%d, Diferença=%f",j,Dif); //Informa
end // Fim do laço de repetição j

mprintf("Bom dia \n Fim da analise..."); // fim

Página 114 – Iterações do sistema através do método Jacobi Gauss-Seidel

//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
//% Resolução de Sistema de Equações Lineares %
//% pelo Método de Jacobe – Gauss Siegel %
//% Dr Alexandre Maniçoba – Labmax.org %
//% Feito para Scilab Versão 2018V100 %
//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clc; // Limpa o console
clear; // Limpa à memória

function [B, w]=iteracao(AA, bb) // Defini a matriz
	[n c] = size (AA);
	for i = 1:n // i variável contadora em n vezes
		bb(i) = bb(i)/AA(i,i);
		AA(i,:) = AA(i,:)/AA(i,i);
	end // Será repetido até aqui
	AA = AA - eye (n,n);
	B = AA; // Atribui o valor AA a variável B
	w = bb; // Atribui o valor bb a variável w
endfunction
	
A = [8 -4 -2; -4 10 -2; -2 -2 10]; // Vetor de Coeficientes
b = [10; 0; 4]; // Vetor de termos independentes
Ep = 10^(-8); // Precisão 1 a 10E-15
it = 100; // Máximo de iterações 1 a 1000
n = length (A(1,:));
x = zeros (n,1); // Define uma matriz de zeros de uma coluna em x
c = ones (n,1); // Define uma matriz de uns de uma coluna em c
k = 0; // Atribui o valor de 0 a variável k
printf("\nA = "); // Imprime na janela do console a mensagem entre aspas
disp (A); // Imprime na janela do console a matriz A
printf("\nb = "); // Imprime na janela do console a mensagem entre aspas
disp(b); // Imprime na janela do console valores de b
printf("\nx = ");// Imprime na janela do console a mensagem entre aspas
disp(x); // Imprime na janela do console valores de x
[B w] = iteracao (A,b); // x = G(x) --> x = Bx+d
printf("\n==================================x(k+1) = w - B.x(k)"); // Imprime na janela do console a mensagem entre aspas
printf("\nw = "); // Imprime na janela do console a mensagem entre aspas
disp(w); // Imprime na janela do console valores de w
printf("\n B = "); // Imprime na janela do console a mensagem entre aspas
disp (B); // Imprime na janela do console valores de B
Bav = abs(spec(B)); // Bav são autovalores de B
[m,n2] = max(Bav); // m = valor absoluto do maior valor de Bav

if m<1 then // Teste de convergência pela norma ||B||<1
	while k<it &abs(max(x-c)) > Ep // Executa até Er<Ep e no máximo it vezes
		c = x; // Atribui os valores de x a variável c
		k = k + 1; // Contador de iterações
		//Processo Iterativo
		printf("\n==================================Iteração %d", k); // Imprime na janela do console a mensagem entre aspas e o valor inteiro de k
		printf("\n-B*w = "); // Imprime na janela do console a mensagem entre aspas
		disp (-B*w); // Imprime na janela do console valores do produto –B*w
		printf("\nx(%d) = ", k); // Imprime na janela do console a mensagem entre aspas com o valor da variável k com números inteiros
		disp(x); // Imprime na janela do console valores de x
		x = w-B*x; //Atribui a x o resultado da equação
		Er(k) = abs(max(x-c)); //Atribui valor absoluto a Er
		printf("\nx(%d) = ", k+1); // Imprime na janela do console a mensagem entre aspas com o valor da variável k+1 com números inteiros
		disp(x); // Imprime na janela do console valores de x
	end // encerra while
	printf("Sistema resolvildo com %d iterações \n\n\n",k);
	// Imprime na janela do console a mensagem entre aspas e valor inteiro de k
	printf("\nA = "); // Imprime na janela do console a mensagem entre aspas
	disp(A); // Imprime na janela do console valores de A
	printf("\nb = "); // Imprime na janela do console a mensagem entre aspas
	disp(b); // Imprime na janela do console valores de b
	printf("\nx = "); // Imprime na janela do console a mensagem entre aspas
	disp(x); // Imprime na janela do console valores de x
	figure (1);
	clf; // Limpa a janela gráfica
	plot(Er); // Plota o gráfico do erro absoluto
	pEr = gca(); // Obtendo o manipulador dos eixos correntes
	pEr.labels_font_size = 5;
	pEr.log_flags = "nln" // Ajustando o eixo X para escala logarítmica
	xtitle("Método de Jacobi e Gauss - Seigel", "Interação", "Erro Relativo (ABS)"); // Imprime na janela do gráfica as mensagem entre aspas: Título, no eixo x e no eixo y
else // Caso as condições não sejam atendidas para convergir
	clc; // Limpa o console
	printf("Não converge para a solução, pois o maior elemento do autovalor de B é maior que 1 \n\n\n"); // Imprime na janela do console a mensagem entre aspas caso a solução não convergir
end // Encerra if/ else

Página 123 – Zero Reais de Funções Reais pelo Método Bissecção

//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
//% Zero Reais de Funções Reais %
//% pelo Método Bissecção %
//% Dr Alexandre Maniçoba – Labmax.org %
//% Feito para Scilab Versão 2018V100 %
//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clf; // Limpa a janela gráfica
clear; // Limpa a memória
clc; // Limpa a janela de console

function [y]=f(x) // Cria uma função f(x)
	y = (-3*x^5)+(2*x^4)+(-5*x^3)+(2*x^2)+x-1 // Função que deseja analisar
endfunction // Fim da função f(x)

function [ym]=grafico(a, b) // Início da função gráfico de (a,b)
	intervalo = b - a; // Atribui o resultado de fiminicial a variável intervalo
	intervalo = intervalo/100; // Reduz o valor do intervalo em 100 vezes
	while(a < b); // Repetição do a enquanto for menor que b
		x = a; // Atribui valor de "a" a x
		y = f(x); // Atribui o valor de f(x) a y
		if y < ym then // Verifica se y é menor que ym
			ym = y; // Se for, atribui o valor de y a ym
		end // Fim da função if
		a = a + intervalo; // Incrementa o valor de "a" com o valor do intervalo
		plot (x, y, '.b'); // Plota um ponto azul nas coordenadas (x, y) no gráfico
	end // Fim do loop while
endfunction // Fim da função gráfico (a,b)

ok = 1; // Atribui o valor 1 a variável ok
while ok == 1 // Laço enqunto ok for igual a 1
	ym = 10000; //atribui 10000 para a variável ym;
	clc; // Limpa a janela de console
	a = input ("Digite o valor de a:"); // Solicita a entrada via teclado e armazena em a
	b = input ("Digite o valor e b:"); // Solicita a entrada via teclado e armazena em b
	ini = a; // Atribui valor de "a" na variável "ini"
	fim = b; // Atribui valor de "b" na variável "fim"
	interm = (a + b)/2; // Atribução da média entre a e b em "interm"
	ym = grafico (a, b); // Evoca a função gráfico() passando a e b e recebendo ym
	// Calcula a raiz da f(x) no intervalo [a, b] com precisão eps1
	x0 = a; // Atribui o valor da variável "a" a variável "x0"
	x1 = b; // Atribui o valor da variável "b" a variável "x1"
	xm = (x0 + x1)./2; // Atribuição da média entre x0 e x1 em "xm" com 2 (./2) algarismos pós vírgula
	eps1 = 10^(-4); // Precisão numérica
	it = 0; // Atribuição do valor 0 a variável it
	
	if f(x0)*f(x1) >= 0 then // Verifica se f(a) e f(b) possuem sinais diferentes
		printf("Valor de f(a) e f(b) devem ter sinal diferente"); // Imprime na janela do console a mensagem entre aspas, caso as funções f(a) e f(b) tiverem sinais iguais
		abort; // Interrompe a avaliação corrente do (if) e retorna
	end; // Fim do if
	
	while abs(f(xm)) > eps1 & it <= 500 // Inicia laço
		if f(x0)*f(xm) > 0 then // Escolhe o intervalo entre "a" e xm ou entre xm e "b"
			x0 = xm; // Atribui o valor de x0 a xm se verdadeiro
		else // Senão
			x1 = xm; // Atribui o valor de x1 xm se falso
		end; // Fim do if
		xm = (x0 + x1)/2; // Atribuição da média entre x0 e x1 em "xm"
		it = it + 1; // Contagem da variável it
		it = it + 1; // Contagem da variável it
	end; // Fim do laço
	
	raiz = xm; // Atribui o valor de xm na variável raiz
	iter = it; // Atribui o valor de it em iter
	if it >= 499 then // Verifica se o número de iterações it é maior que 498
		printf("Não converge nesse intervalo!"); // Informa ao usuário
		abort; // Interrompe o laço ok == 1
	else; // Caso o número de iterações da variável it for menor que 499
		for i = a:((abs(a) + abs(b))/50):raiz // Laço for imprime a linha amarela horizontal do gráfico
		y = 0; // Atribui o valor zero para y
		x = i; // Atribui o valor de i em x
		plot (x, y, '.y'); // Plotagem do gráfico nos eixos x, y em amarelo
		end // Fim do (for)
		for j=0: (ym/30):ym // Laço for imprime a linha amarela vertical do gráfico
		y = j; // Atribui o valor da vaiável contado j a y
		x = raiz; // Atribui o valor da variável raiz a x
		plot(x, y, '.y'); // Plotagem do gráfico nos eixos x, y em amarelo
		end // Fim do (for)
		y = 0; // Atribui o valor zero para x
		x = raiz; // Atribui o valor de raiz na variável x
		plot(x, y , '.r'); // Plotagem em vermelho o ponto formado por x,y
		printf("Raiz Bissecção é %10.15f com %f iterações", raiz, iter); // Imprime na tela da janela do console a mensagem entre aspas com o valor das variáveis raiz e iter em números reais
		ok = input ("1 - Repetir ou [Enter] para sair"); // Imprimi na tela da janela do console a mensagem entre aspas e atribui o valor digitado 1 ou enter na variável ok
		clf; // Limpa a janela gráfica
	end //Fim do if
end // Fim do laço ok == 1

Página 129 – Zero Reais de Funções Reais pelo Método de Newton

//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
//% Zero Reais de Funções Reais %
//% pelo Método de Newton %
//% Dr Alexandre Maniçoba – Labmax.org %
//% Feito para Scilab Versão 2018V100 %
//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clf; //limpa a janela gráfico
clear; //limpa a memoria
clc; //limpa a janela do console
function [y]=f(x) //cria f(x)
	y = (2*x^2) - 4*x - 1; //aqui vai a função
endfunction //fim da função f(x)

function [y]=df(x) // cria uma função diferencial df(x) - Derivada da função f(x)
	y = (4*x) - 4; //aqui vai a derivada de f(x)
endfunction //fim da função df(x)

function [ym]=grafico(a, b) //inicio da função grafico de (a,b)
	intervalo = b - a; //atribui o resultado de fim-ini a variavel intervalo
	intervalo = intervalo/100; //reduz o valor do intervalo em 100 vezes
	while(a < b); //repetição de ini enquanto for menor que b
		x = a; //atribui valor de "a" a x
		y = f(x); //atribui o valor de f(x) a y
		if y < ym then //cerifica se y é menor que ym
			ym = y; //se (for) atrubui o valor de y a ym
		end //fim da função if
		a = a + intervalo; //incrementa o valor de "a" com o valor do intervalo
		plot (x, y, '.b'); //plota um ponto a zul nas coordenadas (x, y) no grafico
	end //fim do loop while
endfunction //fim da função gráfico (a,b)

ok = 1; //atribuia o valor 1 a variável ok
while ok == 1 //laço enquanto ok for igual a 1
	ym = 10000; //atribui 10000 para a variável ym;
	clc; //limpa a janela do console
	a = -1; //atribuição de um valor a variável "a" - "a" é o inicio do intervalo
	b = 1; //atribuição de um valor a variável "b" - "b" é o fim do intervalo
	xi = a; //atribui o valor da variável "a" a variável "xi"
	ini = a; //atribui o valor da variável "a" a variável "ini"
	x(1) = a; //atrubui o valor da variável "a" em x(1)
	ea = 100; //atribui o valor 100 a variável "ea"
	fim = b; //atribui o valor da variável b na varável "fim"
	interm = (a+b)/2; //atribução da media entre a e b em "interm"
	ym = grafico(a,b); //evoca a função gráfico() passando a e b e recebendo ym
	//calcula a raiz da f(x) no intervalo [a, b] com precisão eps1
	x0 = a; //atribui o valor da variável "a" a variável "x0"
	x1 = b; //atribui o valor da variável "b" a variável "x1"
	xm = (x0 + x1)./2; //atribução da media entre x0 e x1 em "xm" com 2 (./2) algarismos pos vírgula
	eps1 = 10^(-6); //precisão
	xf = f(xi); //atribui o valor da função f(xi) em xf
	it = 1; //atribuição do valor 1 a variável it
	if f(x0)*f(x1) >= 0 then //verifica se f(a) e f(b) possuem sinais diferentes
		printf("Valor de f(a) e f(b) devem ter sinal diferente"); //imprime na janela do console a mensagem entre aspas caso as funções f(a) e f(b) tiverem sinais iguais
		abort; //interrompe a avaliação corrente do (if) e retorna
	end; //fim do if

	while abs (ea) >= eps1 & it < 500; //looping enquanto a variável ea de valor absoluto for igual a eps1 e it < 500
		x(it + 1) = x(it) - f(x(it))/df(x(it)); //atribuição da interação de x(it) com a função f(x(it)) e sua derivada df(x(it)), sendo que it a um valor atual, na variável x(it + 1) futura
		ea = abs((x(it+1) - x(it))/x(it+1)*100); //atribui o valor da interação x(it + 1) posterior com x(it) atual na varável ea
		it = it + 1; //contagem da variável it
	end //fim do loop while

	xm = x(it); //atribui o valor x(it) em xm
	raiz = xm; //atribui o valor de xm na variável raiz
	iter = it; //atribui o valor de it em iter
	if it >= 499 then //verifica se o numero de iterações it é maior que 498
		printf("Não converge nesse valor!") //informa o usuário
		abort; //interrompe o laço ok ==1
	else; //caso o número de iterações it for menor que 499
		for i = a:((abs(a) + abs(b))/50):raiz //laço (for) imprime a linha amarela horizontal do gráfico
			y = 0; //atribui o valor zero para y
			x = i; //atribui o valor da variável i em x
			plot (x, y, '.y'); //plotagem do grafico nos eixos x, y em amarelo
		end //fim do for
		for j=0: (ym/30):ym //laço ()for) imprime a linha amarela vertical do gráfico
			y = j; //atribui o valor da vaiável contado j a y
			x = raiz; // atribui o valor da variavel raiz a x
			plot(x, y, '.y'); //plotagem do grafico nos eixos x, y em amarelo
		end //fim do for
		y = 0; //atribui o valor 0 a varável y
		x = raiz; // atribui o valor da variável raiz em x
		plot(x, y , '.r'); //protagem em vermelho no ponte x,y
		printf("Raiz Newton é %10.15f com %f iterações", raiz, iter); //imprime na tela da janela do console a mensagem entre aspas com o valor das variáveis raiz e iter
		ok = input ("1 - Repetir ou [Enter] para sair"); //imprimi na tela da janela do console a mensagem entre aspas e atribui o valor digitado 1 ou enter na variável ok
		clf; //limpa a janela gráfica
	end //fim do if
end // fim do laço ok == 1

Página 142 – Integração numérica – Regra dos Trapézios

//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
//% Integração Numérica %
//% Regra dos Trapézios – Elano Diniz %
//% Dr Alexandre Maniçoba – Labmax.org %
//% Feito para Scilab Versão 2018V100 %
//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear; // Limpa a memoria
clc; // Limpa a janela do console

function [y]=f(x); // Cria a função f(x)
	y = sin(x) + 2; //Especifica a função
endfunction //Encerra a função

a = input ("Entre com a: ") //Pede para o usuário entrar com um valor que será vinculado a variável “a”, em que é valor inicial do intervalo de f(x)
b = input ("Entre com b: ") //Pede para o usuário entrar com um valor que será vinculado a variável “b”, em que é o valor final do intervalo de f(x)
n = input ("Entre com n: ") //Pede para o usuário entrar com um valor que será vinculado a variável “n”, em que n é a quantidade de trapézios a dividir
h = (b-a)/n; //Atribui valor a variável “h” de amplitude
xi = a; //Atribui a f(xi) inicial a variável “a”
soma = 0; //Atribui o valor 0 a variável soma
x = 0; //Atribui o valor 0 a variável soma
c = f(xi); //Atribui f(xi) a variável “c”
xi = a + h; //Atribui ao ponto inicial xi
printf ('\nx%g=%4.3f', x, a); //Imprime na tela do console as variáveis “x” e “a”
printf (' f(x%g)=%4.3f', x, c); //Imprime na tela do console as variáveis “x” e “c”

for i = 1: n-1; //Especifica uma contagem decrescente para repetição do comando
	x = x + 1; //Inicia a contagem da variável x
	c = f(xi); //Atribui f(xi) a variável “c”
	soma = soma + f(xi); //Faz a soma da variável soma junto a função f(xi)
	printf('\nx%g=%4.3f', x, xi); //Imprime na tela do console as variáveis “x” e “xi”
	printf(' f(x%g)=%4.3f', x, c); //Imprime na tela do console as variáveis “x” e “c”
	printf(' soma=%4.3f', soma); //Imprime na tela do console a variável soma
	xi = xi + h; //Faz uma soma que atribui a xi
end //Encerra o comando for

c = f(xi); //Atribui f(xi) a variável “c”
I = (h/2)*(f(a) + f(b) + 2*soma); //Faz a soma de todas as sub interações “n”, resultando na área
printf('\nx%g=%4.3f', n, b); //Imprime na tela do console as variáveis “n” e “b”
printf(' f(x%g)=%4.3f', n, c); //Imprime na tela do console as variáveis “n” e “c”
printf(' soma=%4.3f', soma); //Imprime na tela do console a variável soma
printf("\n\n a=%4.3f b=%4.3f n=%4.3f", a, b, n); //Imprime na tela do console as variáveis “a”, “b” e “n”
printf("\n(b-a)=%4.3f", (b-a)); //Imprime na tela do console o intervalo entre “a” e “b”
printf(' h=%4.3f', h); //Imprime na tela do console a variável “h”
printf(" I=%4.3f", I); //Imprime na tela do console a variável “I” que equivale a área delimitada no intervalo[a,b] da função f(x)

Página 145 – Integração numérica – Regra de Simpson

//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
//% Integração Numérica %
//% Regra de Simpson %
//% Dr Alexandre Maniçoba – Labmax.org %
//% Feito para Scilab Versão 2018V100 %
//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear; //Limpa memória
clc; //Limpa a janela do console

function [y]=f(x)//Cria a função f(x)
	y = sin(x)+2 //Define a função f(x)
endfunction //Encerra a função

a = input ("Entre com a: ") //Imprime na tela e atribui valor digitado pelo usuário a variável a (intervalo)
b = input ("Entre com b: ") //Imprime na tela e atribui valor digitado pelo usuário a variável b(intervalo)
n = input ("Entre com n: ") //Imprime na tela e atribui valor digitado pelo usuário a variável n(subintervalos)
h = (b-a)/(2*n) //Atribui valor da amplitude a variável "h"
soma = f(a) + f(b) //Soma das duas função no intervalo [a,b]
for i=1 : n-1 //Inicia um contagem com decremento -1
	soma=soma + 2*f(a+h*2*i) + 4*f(a+h*(2*i-1)) //Realiza a soma de todas as subdivisões de n e atribui a variável soma
end //Encerra a função for
soma = soma + 4*f(a+(2*n-1)*h) //atribui valor a variável soma
I = (h/3)*soma //Atribui o valor da soma da área final a I
printf("\n a = %4.3f b = %4.3f n = %4.3f",a,b,n) //imprime na tela as várias a,b,n
printf("\n (b-a) = %4.3f", (b-a))//imprime na tela o intervalo entre a e b
printf(' h = %4.3f',h)//imprime na tela amplitude h
printf(" I = %4.3f",I) //imprime na tela o valor da variável I da área.

Página 154 – Método dos Quadrado I – Ajuste de Curvas

//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
//% Metodo dos Quadrado I - Ajuste de Curvas %
//% Ajuste de Curva de Tendencia - D. Justo %
//% Dr Alexandre Maniçoba – Labmax.org %
//% Referência Dagoberto Justo %
//% https://www.youtube.com/watch?v=Tm9wH5dUBIY %
//% Feito para Scilab Versão 2018V100 %
//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear; //Limpa a memória
clc; //Limpa a janela de console
x = [1 2 3 4 5 6 7 8 9 10]; //Pontos da coordenada x
y = [1 1 2 4 5 6 3 8 9 7]; //Pontos da coordenada y
n = size(x,1); //Ajusta a quantidade de números de pontos (tamanho do vetor)
M = [n sum(x); sum(x) sum(x.^2)] //Faz a soma das linha e colunas da matriz M
w = [sum(y); sum(x.*y)] //Faz a soma das linha e colunas da matriz w
a = inv(M)*w; //Atribui a matriz inversa M*w a variável "a""
X = 0:0.1:n; //Atribui valores ao vetor "x"
Y = a(1) + a(2)*X //Atribui os valores do vetor "b" em função de f(x)
clf; //Limpa a janela gráfica
plot(x,y,'b*'); //Plota os pontos (x,y) em azul
ps = input('Pressione Enter'); //Aguarda até o pressionamento de enter
plot(x,y, 'c');//Plota a curva em azul claro
ps = input('Pressione Enter'); //Aguarda até o pressionamento de enter
plot(X,Y,'r');//Plota a curva tendência em vermelho

Página 157 – Metodo dos Quadrado II – Ajuste de Curvas

//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
//% Metodo dos Quadrado II - Ajuste de Curvas %
//% Ajuste de Curva de Tendencia - D. Justo %
//% Dr Alexandre Maniçoba – Labmax.org %
//% Referência Dagoberto Justo %
//% https://www.youtube.com/watch?v=Tm9wH5dUBIY %
//% Feito para Scilab Versão 2018V100 %
//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear; //Limpa a memória
clc; //Limpa a janela de console
x = [1 2 3 4 5 6 7 8 9 10]; //Pontos da coordenada x
y = [1 1 2 4 5 6 3 8 9 7]; //Pontos da coordenada y
n = size(x,1); //Ajusta a quantidade de números de pontos(tamanho do vetor)
clf; // Limpa a janela gráfica
plot(x,y, 'b*'); // Plota os pontos (x,y) em azul
ps = input('Pressione Enter'); // Aguarda até pressionamento de enter
plot(x,y, 'c'); // Plota a curva em azul claro
pmax = 5; //atribui valor 5 para pmax, que é o grau do polinomio maximo
for p=1:pmax //Faz uma contagem de 1 a 5, que é o valor de pmax do grau polinomial
	for i=1:p+1 //Faz uma contagem implementando +1 i linhas
		for j=1:p+1 //Faz uma contagem implementando +1 j colunas
			M(i,j)=sum(x.^(i+j-2)); //Faz a soma das linhas e colunas da matriz M
		end // Encerra o comando for
	end // Encerra o comando for
	for i=1:p+1 //Faz uma contagem de implemento +1 das linhas
		w(i) = sum(y.*x.^(i-1)); //soma da linha do vetor w
	end // Encerra o comando for
	a = inv(M)*w; //Calcula a matriz inversa M*w e atribui a variável “a”
	X = 0:0.1:10; //Atribui os valores de 0 a n com 0.1 ao vetor X
	Y = 0; //Atribui valor 0 a função de f(x)
	R = 0; // Atribui valor 0 a variável residual
	for i=1:p+1 //faz a contagem do grau polinomial f(x)
		Y = Y+a(i)*X.^(i-1); //atribui o valor da equação polinomial a f(x)
		R = R+a(i)*x.^(i-1); //atribui valor ao resíduo
	end //Encerra o comando for
	R = R - y; //atribui o valor de decremento do resíduo anterior com o ponto do gráfico
	clf; //Limpa a janela gráfica
	plot(x,y,'b*'); //Plota os pontos (x,y) em azul
	plot(x,y,'c'); // Plota a curva em azul claro
	plot(X,Y,'r'); // Plota a curva de tendência em vermelho
	printf('f(x) grau %d R = %f4.2',p,R); //imprime na janela do console o numero do grau "p" e o valor do resíduo "R"
	ps = input (''); //Aguarda até pressionamento de enter
end //Encerra o comando for

Página 169 – Algoritmo de Lógica Paraconsistente. Estado Lógico Paraconsistente

//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
//% Algoritmo de Lógica Paraconsistente %
//% Estado Lógico Paraconsistente %
//% Autor: Arnaldo de Carvalho Junior %
//% Junho 2020 %
//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
// Enquanto resposta for 's', realiza análise LPA2v
resp='s'
while resp=='s'||resp=='S'
// Limpa memória
clear
clc

// Exibe mensagens e aguarda entrada de dados na janela de console
disp "*** LÓGICA PARACONSISTENTE ANOTADA COM ANOTAÇÃO DE 2 VALORES"
disp " LPA2v"
disp ""
Mi = input("Informe o Grau de Evidência Favorável (entre 0.00 e 1.00): ")
disp ""
Lbda=input("Informe o Grau de Evidência Desfavorável (entre 0.00 e 1.00): ")
disp ""
Lcc=input("Informe o Controle de Certeza ({0,1} Valor default = 0.5}")
disp ""
Lcct=input("Informe o Controle de Contradição ({0,1} Valor default = 0.5}")
clc
disp "*** LÓGICA PARACONSISTENTE ANOTADA COM ANOTAÇÃO DE 2 VALORES"
disp ""
disp " LPA2v"
disp ""
mprintf("Grau de Evidência Favorável (Mi) = %7.2f \n",Mi)
mprintf("Grau de Evidência Desfavorável (Lambda) %7.2f \n",Lbda)
disp ""

// Executa análise LPA2v
GC=Mi-Lbda; // Cálculo GC
mprintf("Grau de Certeza (GC) = %7.2f \n",GC);
GCT=Mi+Lbda-1; // Cálculo GCT
mprintf("Grau de Contradição (GCT) = %7.2f \n",GCT);
disp ""
D=((1-abs(GC))^2+GCT^2)^0.5; // Cálculo Distância

if D>1 then D=1;
else D=D;
end

if GC>0 then GCR=1-D;
else if GC<0 then GCR=D-1;
else GCR=0.5;
end
end
MIE=(GC+1)/2; // Cálçulo MIE
MIER=(GCR+1)/2; // Cálculo MIER
mprintf("Grau de Evidência Resultante (MIE) = %7.2f \n",MIE);
mprintf("Grau de Evidência Resultante Real (MIER) = %7.2f \n",MIER);
MIEct=(GCT+1)/2; // Cálculo MIEct
mprintf("Grau de Contradição Resultante (MIEct) = %7.2f \n",MIEct);
disp ""

// Classificação dos 12 estados lógicos paraconsistentes
if GC>=Lcc then S1="Verdadeiro (V)"
else if GC<=-Lcc then S1="Falso (F)"
else if GCT>=Lcct then S1="Inconsistente (T)"
else if GCT<=-Lcct then S1="Indeterminado (I)"
else if GC>=0 && GC<Lcc && GCT>=0 && GCT<Lcct && GC>=GCT
then S1="Quase Verdadeiro tendendo a Inconsistente (Qv=>T)"
else if GC>=0 && GC<Lcc && GCT>=0 && GCT<Lcct && GC<GCT
then S1="Inconsistente tendendo a Verdadeiro (T=>V)"
else if GC>=0 && GC<Lcc && GCT>-Lcct && GCT<=0 &&
GC>=abs(GCT) then S1="Quase Verdadeiro tendendo a Indeterminado (Qv=>I)"
else if GC>=0 && GC<Lcc && GCT>-Lcct && GCT<=0 && GC<abs(GCT) then S1="Indeterminado tendendo a Verdadeiro (I=>V)"
else if GC>-Lcc && GC<=0 && GCT>-Lcct && GCT<=0 && abs(GC)>=abs(GCT) then S1="Quase Falso tendendo a Indeterminado (Qf=>I)"
else if GC>-Lcc && GC<=0 && GCT>-Lcct && GCT<=0 && abs(GC)<abs(GCT) then S1="Indeterminado tendendo a Falso (I=>F)"
else if GC>-Lcc && GC<=0 && GCT>=0 && GCT<Lcct && abs(GC)>=GCT then S1="Quase Falso tendendo a Inconsistente (Qf=>T)"
else S1="T=>F"
end
end
end
end
end
end
end
end
end
end
end
// Apresenta Resultados LPA2v
mprintf("Estado lógico Paraconsistente é:");
disp(S1)

//getf(LPA2v_Lattice_scilab.sce)
//-------------------------------------------------------
// MONTA GRÁFICO
// ------------------------------------------------------
// desenha o reticulado
clf
gce().font_size = 2;
plot([-1 0],[0 1],'k','LineWidth',2)
plot([0 1], [1 0],'k','LineWidth',2)
plot([1 0], [0 -1],'k','LineWidth',2)
plot([0 -1],[-1 0],'k','LineWidth',2)
plot([-1 1],[0 0],'k','LineWidth',2)
plot([0 0], [-1 1],'k','LineWidth',2)
// Desenha Controles de Certeza e de Contradição
plot([Lcc,Lcc],[(1-Lcc),(-1+Lcc)],'r','LineWidth',2)
plot([-Lcc,-Lcc],[(1-Lcc),(-1+Lcc)],'r','LineWidth',2)
plot([-Lcc,Lcc], [(1-Lcc),(1-Lcc)],'r','LineWidth',2)
plot([-Lcc,Lcc],[(-1+Lcc),(-1+Lcc)],'r','LineWidth',2)
plot([-Lcc,Lcc],[(-1+Lcc),(1-Lcc)],'r','LineWidth',2)
plot([-Lcc,Lcc],[(1-Lcc),(-1+Lcc)],'r','LineWidth',2)

//%%%%%%%%%%%%%%% inserindo valores no reticulado%%%%%%%%%%%%%%
//-------------------------------------------------------
// Valores relacionados a Lambda
// Monta a reta vertical para o valor de Lambda até encostar no reticulado
//% line([Lambda Lambda],[-1.1 (-1-Lambda)],'Color','r','LineWidth',1,'LineStyle',':');
plot([-Lbda -Lbda],[-1 (-1+Lbda)],'g','LineWidth',2,'LineStyle',':')

// monta a reta diagonal dentro do reticulado
plot([-Lbda (1-Lbda)],[(-1+Lbda)
Lbda],'g','LineWidth',2,'LineStyle',':')

// Valores relacionados a Mi
// monta a reta vertical para o valor de Mi até encostar no reticulado
plot([Mi Mi],[-1 (-1+Mi)],'b','LineWidth',2,'LineStyle',':')

// monta a reta diagonal dentro do reticulado
plot([Mi (-1+Mi)],[(-1+Mi)Mi],'b','LineWidth',2,'LineStyle',':')
// monta a reta diagonal dentro do reticulado
plot([Mi (-1+Mi)],[(-1+Mi)Mi],'b','LineWidth',2,'LineStyle',':')
//-------------------------------------------------------
// exibe valores
xstring(Mi,-0.95,['<=','$\mu$'])
xstring(-Lbda,-0.95,["<=",'$\lambda$'])
xstring(GC, (GCT-0.06),["<=",'$\epsilon$','(',string(GC),',',string(GCT),")"])
xstring(0.3,0.9,"GRAU DE EVIDENCIA ")
xstring(0.3,0.8,["RES. REAL = ",string(MIER)])
xstring(GCR,-0.06,"|")
xstring((GCR-0.03),-0.12,"GCR")
xtitle( 'Diagrama Lattice LPA2v', boxed = %t)
disp(' ')
resp=input('Deseja realizar novo estudo (S/N) ?',"s")
R=isempty(resp)
if R==1 || resp<>"s" then resp = 'N'
else resp='S'
end
disp(' ')
disp('ATÉ BREVE!');
end;

Página 182 – Algoritmo de Indicador de Qualidade do Canal

//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
//% Algoritmo de Indicador de Qualidade do Canal %
//% CQIpal2v %
//% Autor: Arnaldo de Carvalho Junior %
//% Junho 2020 %
//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
// Se digitar "s" executa uma nova análise do CQIpal2v
resp='s'
while resp=='s'||resp=='S'
// Limpa tela e variáveis
clear
clc
// Monta tela inicial e aguarda entrada de valores
disp "*** INDICADOR DA QUALIDADE DO CANAL DE TV DIGITAL
LPA2v ***"
disp ""
PMi= input("Informe a Potência do Canal (0 e 100%): ")
MI=PMi/100;
disp ""
QLbda=input("Qual o Índice de Qualidade do Canal (0 e 100%): ")
LBDA=1-QLbda/100;
// Limpa a tela e apresenta a análise LPA2v
clc
disp "*** INDICADOR DA QUALIDADE DO CANAL DE TV DIGITAL
LPA2v"
disp ""
mprintf("Potência do Canal = %7.2f \n",PMi,'por cento')
mprintf("Qualidade do Canal = %7.2f \n",QLbda,'por cento')
mprintf("Grau de Evidência Favorável (Mi) = %7.2f \n",MI)
mprintf("Grau de Evidência Desfavorável (Lambda) %7.2f\n",LBDA)
disp ""
GC=MI-LBDA; // Cálculo GC
mprintf("Grau de Certeza (GC) = %7.2f \n",GC);
GCT=MI+LBDA-1; // Cálculo GCT
mprintf("Grau de Contradição (GCT) = %7.2f \n",GCT);
disp ""
D=((1-abs(GC))^2+GCT^2)^0.5; // Cálculo Distância
if D>1 then D=1;
else D=D;
end
if GC>0 then GCR=1-D;
else if GC<0 then GCR=D-1;
else GCR=0.5;
end
end
MIE=(GC+1)/2; // Calculo MIE
mprintf("Grau de Evidência Resultante (MIE) = %7.2f
\n",MIE);
MIct=(GCT+1)/2; // Calculo MIEct
mprintf("Grau de Contradição Resultante (MIct) = %7.2f
\n",MIct);
MIER=(GCR+1)/2; // Calculo MIER
mprintf("Grau de Evidência Resultante Real (MIER) = %7.2f
\n",MIER);
if MIER>0.9 then S1="EXCELENTE"
else if MIER>0.7 && MIER<=0.9 then S1="BOM"
else if MIER>0.5 && MIER<=0.7 then S1="REGULAR"
else if MIER>0.3 && MIER<=0.5 then S1="FRACO"
else if MIER>0.1 && MIER<=0.3 then
S1="PÉSSIMO"
else S1="SEM SINAL"
end
end
end
end
end
// Apresenta resultado do CQIpal2v
disp ""
mprintf("O Indicador da Qualidade Final do Canal é :")
disp(S1)
disp("")
resp=input("Deseja realizar novo estudo (S/N) ?","s")
R=isempty(resp)
if R==1 || resp<>"s" then resp = 'N'
else resp="S"
end
disp("")
end
disp('ATÉ BREVE!');

Página 196 – Algoritmo de Link Quality Estimator para IWSN

//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
//% Algoritmo de Link Quality Estimator para IWSN %
//% Autor: Arnaldo de Carvalho Junior %
//% Junho 2020 %
//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
// Limpa memória
clear;
clc;
// Carrega arquivo de Medições
M=csvRead('C:\LabMax\medidas.txt',[],[],'string'); // Escrever o caminho inteiro do diretório correto para o arquivo de medidas.
// Converte vetores de string em vetores numéricos LQI e RSSI
LQI = part([M], 7:9);
LQI = strtod(LQI)
RSSI = part([M],19:21);
RSSI = strtod(RSSI);

// Aplica normalização de RSSI e LQI
L=length(RSSI);
for x=1:1:L;
if RSSI(x,1)>90 then
MIrssi(x,1)=1;
elseif RSSI(x,1)<10 then
MIrssi(x,1)=0;
else MIrssi(x,1)=(RSSI(x,1)-10)/80; // Linear 1 >= 90% e 0 <= 10%
end
MIlqi(x,1)=0.404*LQI(x,1)+60.134-(227/LQI(x,1))-19*log(LQI(x,1)); // Conversão Não Linear de LQI
if MIlqi(x,1)>1 then
MIlqi(x,1)=1;
elseif MIlqi(x,1)<0 then
MIlqi(x,1)=0;
else
end
end

// Rede Neural Paraconsistente para LQEpal2v
// Filtro LPA2v LQI
FL1=0.25; // Define Valor do Fator de Aprendizagem do Filtro LPA2v para RSSI
FL2=0.25; // Define Valor do Fator de Aprendizagem do Filtro LPA2v para LQI
MI1=MIrssi;
MI2=MIlqi;
MIE11(1,1)=MI1(1,1); MIE12(1,1)=MIE11(1,1);
MIE13(1,1)=MIE12(1,1); // Valores Iniciais CNAPapctx Filtro 1
MIE21(1,1)=MI2(1,1); MIE22(1,1)=MIE21(1,1);
MIE23(1,1)=MIE22(1,1); // Valores Iniciais CNAPapctx Filtro 2
for x=2:1:L;
// Filtro LPA2v RSSI
MIE11(x,1)=MI1(x,1)*FL1+MIE11(x-1,1)*(1-FL1);
MIE12(x,1)=MIE11(x,1)*FL1+MIE12(x-1,1)*(1-FL1);
MIE13(x,1)=MIE12(x,1)*FL1+MIE13(x-1,1)*(1-FL1);
// Filtro LPA2v LQI
MIE21(x,1)=MI2(x,1)*FL2+MIE21(x-1,1)*(1-FL2);
MIE22(x,1)=MIE21(x,1)*FL2+MIE22(x-1,1)*(1-FL2);
MIE23(x,1)=MIE22(x,1)*FL2+MIE23(x-1,1)*(1-FL2);
// NAP de Saída LQEpal2v
MI(x,1)=MIE13(x,1);
LBDA(x,1)=1-MIE23(x,1);
GC(x,1)=MI(x,1)-LBDA(x,1);
GCT(x,1)=MI(x,1)+LBDA(x,1)-1;
D(x,1)=((1-abs(GC(x,1)))^2+GCT(x,1)^2)^0.5;
if D(x,1)>1 then
D(x,1)=1;
else
end
if GC(x,1)>0 then
GCR(x,1)=1-D(x,1);
elseif GC(x,1)<0 then
GCR(x,1)=D(x,1)-1;
else
GCR(x,1)=0;
end
MIER(x,1)=(GCR(x,1)+1)/2;
end
LQEpal2v=MIER;
// Gera Gráficos Saída
// Gráficos RSSI
subplot(221)
plot(MIrssi,'r','LineWidth',2); // plots MIrssi
xgrid(1);
title("RSSI Normalizado",'fontsize',3)
subplot(223)
plot(MIE13,'r','LineWidth',2); // plots MIrssi Filtrado
xgrid(1);
title("RSSI Normalizado e Filtrado",'fontsize',3)

// Gráficos LQI
show_window(1)
subplot(221)
plot(MIlqi,'g','LineWidth',2);// plots MIlqi
xgrid(1);
title("LQI Normalizado",'fontsize',3)
subplot(223)
plot(MIE23,'g','LineWidth',2)// plots MIlqi Filtrado
xgrid(1);
title("LQI Normalizado e Filtrado",'fontsize',3)

// Gráfico LQEpal2v
show_window(2)
plot([0,L],[1,1],'k','LineWidth',1,'LineStyle','-');
plot([0,L],[0.9,0.9],'k','LineWidth',2,'LineStyle','-');
plot([0,L],[0.7,0.7],'k','LineWidth',2,'LineStyle','-');
plot([0,L],[0.5,0.5],'k','LineWidth',2,'LineStyle','-');
plot([0,L],[0.3,0.3],'k','LineWidth',2,'LineStyle','-');
plot([0,L],[0.1,0.1],'k','LineWidth',2,'LineStyle','-');
plot(LQEpal2v,'b','LineWidth',2);// plots LQEpal2v
xstring(0.3,0.9,"EXCELENTE");
xstring(0.3,0.7,"BOM");
xstring(0.3,0.5,"REGULAR");
xstring(0.3,0.3,"FRACO");
xstring(0.3,0.1,"PÉSSIMO");
xstring(0.3,0.0,"SEM CONEXÃO");
title("LQEpal2v",'fontsize',3)
xgrid(1);

Patch Antenna Calculator – Bonus

//********************************************************************
//      Patch Antenna Calculator - PAC v2.00                         *
//      Written by : Dr. Alexandre Manicoba De Oliveira &            *
//                   LabMax, IFSP - Brazil                           *
//       Based in Matlab Program Matlab of the Sung-Woo Lee &        *
//                    Zhiyong Huang, Arizona State University        *
//     Aug. 19, 2023                                                 *
//********************************************************************

// Subfunctions
// ************
function y=sincc(x)
  // normalized sinc function, sin(pi*x)/(pi*x), no checks on the input
  y = sin(pi*x)./(pi*x);
  y(x==0) = 1;
endfunction

function [Ethval, Eth]=E_th(phir, h, ko, Leff, Emax)
  ARG=cos(phir).*h.*ko./2;
  Ethval=(sincc(ARG./pi).*cos(sin(phir).*ko*Leff./2))./Emax;
  Eth=20*log10(abs(Ethval));
  Eth(phir>pi/2&phir<3*pi/2)=-60;
  Eth(Eth<=-60)=-60;
endfunction

function [Ephval, Eph1]=E_ph(thr, h, ko, W, Emax)
  ARG1=sin(thr).*h.*ko./2;
  ARG2=cos(thr).*W.*ko./2;
  Ephval=sin(thr).*sincc(ARG1./pi).*sincc(ARG2./pi)./Emax;
  Eph1=20.0*log10(abs(Ephval));
  Eph1(Eph1<=-60)=-60;
endfunction

//function [D,DdB]=dir_rect(W,h,Leff,~,ko)
function [D, DdB]=dir_rect(W, h, Leff, ko)  
  th=0:180; phi=[0:90 270:360];
  [t,p]=meshgrid(th.*pi/180,phi.*pi/180);
  X=ko*h/2*sin(t).*cos(p);
  Z=ko*W/2*cos(t);
  Et=sin(t).*sincc(X/pi).*sincc(Z/pi).*cos(ko*Leff/2*sin(t).*sin(p));
  U=Et.^2;
  dt=(th(2)-th(1))*pi/180;
  dp=(phi(2)-phi(1))*pi/180;
  Prad=sum(sum(U.*sin(t)))*dt*dp;
  D=4.*pi.*max(max(U))./Prad;
  DdB=10.*log10(D);
endfunction

function [G1, G12]=sintegr(W, L, ko)
  th=0:1:180; t=th.*pi/180;
  ARG=cos(t).*(ko*W/2);
  res1=sum(sincc(ARG./pi).^2.*sin(t).^2.*sin(t).*((pi/180)*(ko*W/2)^2));
  res12=sum(sincc(ARG./pi).^2.*sin(t).^2.*besselj(0,sin(t).*(ko*L)).*sin(t).*((pi/180)*(ko*W/2)^2));
  G1=res1./(120*pi^2); G12=res12./(120*pi^2);
endfunction

//////////////////////////////////////
// Subfunction
//function [D,DdB]=dir_cir(~,ae,ko)
function [D, DdB]=dir_cir(ae, ko)
  th=0:90; phi=0:360;
  [t,p]=meshgrid(th.*pi/180,phi.*pi/180);
  x=sin(t).*ko.*ae;
  J0=besselj(0,x); J2=besselj(2,x);
  J02P=J0-J2; J02=J0+J2;
  Ucirc=(J02P.*cos(p)).^2 + (J02.*cos(t).*sin(p)).^2;
  Umax=max(max(Ucirc));
  Ua=Ucirc.*sin(t).*(pi./180).^2;
  Prad=sum(sum(Ua));
  D=4.*pi.*Umax./Prad;
  DdB=10.*log10(D);
endfunction

clear all;
warning off;
option=1;
filename=[];
patchm=1;
warning on;
pi=%pi;
    // Input Parameters (freq, epsr, height, Yo)
    freq=[];
    // 0<freq<10.1
    while (isempty(freq)|freq>10.1)
       clc;
       disp('LabMax - Patch Antenna ToolBox V. 2.00   ');
       freq=input('Enter the Frequency (GHz) = ');
       //freq=1.5; 
    end

    er=[];
    while (isempty(er)|er>10.1)
       // 0<er<10.1
       er=input('Enter the Dielectric Constant of Substrate = ');
       //er=4.3;
    end

    h=[];
    while (isempty(h)|h>5.1)
       // 0<h<5.1
       h=input('Enter the Height of the Substrate (mm) = ');
       //h=1.6;
       h=h/10;
    end
    option1=2;
    if option1==1
        Yo=[];
        while isempty(Yo)
            Yo=input(['input the position of the recessed feed point ' ...
                      'relative to the leading radiating edge (mm) = ']);
            Yo=Yo/10;
        end
    else
        Zin=[];
        while (isempty(Zin)|Zin>75.1)
            // 0<Zin<75.1
            Zin=input(['Input the desired input impedance Zin (ohms) = ']);
            //Zin=50;
        end
    end

  // Compute W, ereff, Leff, L (in cm)
  W=30.0/(2.0*freq)*sqrt(2.0/(er+1.0));
  ereff=(er+1.0)/2.0+(er-1)/(2.0*sqrt(1.0+12.0*h/W));
  dl=0.412*h*((ereff+0.3)*(W/h+0.264))/((ereff-0.258)*(W/h+0.8));
  lambda_o=30.0/freq;
  lambda=30.0/(freq*sqrt(ereff));
  Leff=30.0/(2.0*freq*sqrt(ereff));
  L=Leff-2.0*dl;
  ko=2.0*pi/lambda_o;
  Emax=sincc(h*ko/2.0/pi);

  //Compute Wfeed of Microstrip Line for Zin impedance 
  Wfeed=((7.48*h*10)/(exp(Zin*(sqrt(er+1.141)/87))))-1.25*20e-6;

  // Normalized radiated field
  //         E-plane pattern : 0 < phi < 90    ;    270 < phi < 360
  //         H-plane pattern : 0 < th < 180
  phi=0:360; phir=phi.*pi./180; Eth=E_th(phir,h,ko,Leff,Emax);
  th=0:180; thr=th.*pi/180.0;
  Eph1=E_ph(thr,h,ko,W,Emax);
  Eph(1:91)=Eph1(91:181); Eph(91:270)=Eph1(181); Eph(271:361)=Eph1(1:91);
  
  // Figure 2
  // ********
  figure(2);
  polarplot(((phi)*pi/180),Eth,style=5);
  title('Radiation Pattern (Theta)','fontsize',[4]);

  figure(3);
  polarplot(((phi+90)*pi/180),Eph,style=5);
  title('Radiation Pattern (Phi)','fontsize',[4]);

  // E-plane HPBW and  H-plane HPBW
  // ******************************
  an=phi(Eth>-3);
  an(an>90)=[];

  EHPBW=2*abs(max(an));
  HHPBW=2*abs(90-min(th(Eph1>-3)));

  // Directivity
  [D,DdB]=dir_rect(W,h,Leff,ko);

  // Input Impedance at Y=0 and Y=Yo
  [G1,G12]=sintegr(W,L,ko);
  Rin0P=(2.*(G1+G12))^-1;
  Rin0M=(2.*(G1-G12))^-1;
  if option1==1
      RinYoP=Rin0P*cos(pi*Yo/L)^2;
      RinYoM=Rin0M*cos(pi*Yo/L)^2;
  else
      YP=acos(sqrt(Zin/Rin0P))*L/pi;
      YM=acos(sqrt(Zin/Rin0M))*L/pi;
  end

  // Display (rectangular)
  option1=2;
  //if(option_a==2)
  //   diary(filename);
  //end
  disp('INPUT PARAMETERS ============================================='); 
  mprintf('Resonant frequency = %4.2f GHz\n',freq); 
  mprintf('Dielectric constant of the substrate = %4.2f\n',er);
  mprintf('Height of the substrate = %4.2f mm\n',(h*10));
  if option1==1
      mprintf('Position of the recessed feed point  = %4.2f mm\n',Yo*10);
  else
      mprintf('Desired resonant input impedance = %4.2f ohms\n\n', Zin);
  end
  disp('OUTPUT PARAMETERS ============================================');
  
  Wa=W*10;      //Patch width in mm
  Waf=Wfeed;    //Microstrip line width in mm
  La=L*10;      //Patch length in mm
  Aa=Wa*2;      //Substrate width and length in mm
  Laf=L*2;      //Slot feed length in mm
                //Slot feed width is 2 mm
                
  mprintf('Microstrip line feed width (Wfeed) = %4.2f mm\n',(Wfeed));
  mprintf('Physical width (W) of patch = %4.2f mm\n',(W*10));
  mprintf('Effective length (Leff) of patch = %4.2f mm\n',(Leff*10));
  mprintf('Physical length (L) of patch = %4.2f mm\n',(L*10));
  mprintf('E-PLANE HPBW = %4.2f degrees\n',EHPBW);
  mprintf('H-PLANE HPBW = %4.2f degrees\n',HHPBW);
  mprintf('Antenna ABS directivity = %4.2f\n',D);
  mprintf('Antenna directivity = %4.2f dB\n',DdB);
  mprintf('\n');
  mprintf('                            %3.1fmm                   \n', Aa);
  mprintf('        +---------------------------------------------+\n');
  mprintf('        |                                             |\n');
  mprintf('        |                                             |\n');
  mprintf('        |                                             |\n');
  mprintf('        |                   %3.1fmm                    |\n',(W*10));
  mprintf('        |         +-----------------------+           |\n');
  mprintf('        |         |                       |           |\n');
  mprintf('        |         |                       |           |\n');
  mprintf('        |         |                       |           |\n');
  mprintf(' %4.1fmm|         |                       |%3.1fmm     |\n',Aa, (L*10));
  mprintf('        |         |      2mm     2mm      |           |\n');
  mprintf('        |         |      +-+     +-+      |           |\n');
  mprintf('        |         |      | |     | |%3.1fmm |           |\n',Laf);
  mprintf('        |         +------+ |     | +------+           |\n');
  mprintf('        |                  |     |                    |\n');
  mprintf('        |                  |%3.1fmm|                    |\n',Wfeed);
  mprintf('        |                  |     |                    |\n');
  mprintf('        +------------------+-----+--------------------+\n');
  
  // Figure 3 - Design of Patch Antenna
  // ********
  figure(4);

  //                        Aa              
  // +------------------------------------------------+
  // |                                                | 
  // |                                                | 
  // |                      Wa                        | 
  // |           A------------------------B           | 
  // |           |                        |           |
  // |           |                        |           |
  // |           |                        |           |
  // |           |                        | La        | Aa
  // |           |                        |           |
  // |           |         2     2        |           |
  // |           |        J-M   N-L       |           |
  // |           |        | |   | | Laf   |           |
  // |           D--------I-F---G-K-------C           |
  // |                      |   |                     |
  // |                      |   |                     |
  // |                      |Waf|                     |
  // +----------------------E---H---------------------+

  //Xa(0)=Aa/2;          Ya(0)=Wa; <-- coord. to center

  Xa(1)=Wa/2;          Ya(1)=(Aa-La)/2+La;  //point A
  Xa(2)=Wa/2+Wa;       Ya(2)=Ya(1);         //point B
  Xa(3)=Xa(2);         Ya(3)=(Aa-La)/2;     //point C
  Xa(4)=Xa(1);         Ya(4)=Ya(3);         //point D

  Xa(5)=Wa-(Waf/2);    Ya(5)=0;             //point E
  Xa(6)=Xa(5);         Ya(6)=(Aa-La)/2;     //point F
  Xa(7)=Wa+(Waf/2);    Ya(7)=Ya(6);         //point G 
  Xa(8)=Xa(7);         Ya(8)=0;             //point H

  Xa(9)=Wa-(Waf/2)-2;  Ya(9)=Ya(6);         //point I
  Xa(10)=Xa(9);        Ya(10)=(Aa-La)/2+Laf;//point J 

  Xa(11)=Wa+(Waf/2)+2; Ya(11)=Ya(6);        //point K
  Xa(12)=Xa(11);       Ya(12)=Ya(10);       //point L

  Xa(13)=Xa(6);        Ya(13)=Ya(10);       //point M
  Xa(14)=Xa(7);        Ya(14)=Ya(10);       //point N

  //               LS                   LS
  //   +------------------------+     +----+
  //   |                        |     |    |
  //   |                        |    L| 2  |L
  //   |                        |    E|  3 |D 
  // L |            1           | L   |    |
  // E |                        | D   +----+
  //   |                        |       LI
  //   |        +-+   +-+       |
  //   |        |2|   |3|       |
  //   +--------+a+-b-+c+-------+
  //      LI      |   | 
  //              | 4 |  
  //              |   |
  //              +---+

  // 1st Retangle - Patch
  // LS
  plot([Xa(1) Xa(2)], [Ya(1) Ya(2)],'b','LineWidth',1)
  // LE
  plot([Xa(1) Xa(4)], [Ya(1) Ya(4)],'b','LineWidth',1)
  // LD
  plot([Xa(2) Xa(3)], [Ya(2) Ya(3)],'b','LineWidth',1)
  // LI
  plot([Xa(4) Xa(3)], [Ya(4) Ya(3)],'b','LineWidth',1)

  // 2nd Retangle - Feed slot
  // LS
  plot([Xa(10) Xa(13)], [Ya(10) Ya(13)],'b','LineWidth',1)
  // LE
  plot([Xa(10) Xa(9)], [Ya(10) Ya(9)],'b','LineWidth',1)
  // LD
  plot([Xa(13) Xa(6)], [Ya(13) Ya(6)],'b','LineWidth',1)
  // LI
  plot([Xa(9) Xa(6)], [Ya(9) Ya(6)],'b','LineWidth',1)

  // 3rd Retangle - Feed slot
  // LS
  plot([Xa(14) Xa(12)], [Ya(14) Ya(12)],'b','LineWidth',1)
  // LE
  plot([Xa(14) Xa(7)], [Ya(14) Ya(7)],'b','LineWidth',1)
  // LD
  plot([Xa(12) Xa(11)], [Ya(12) Ya(11)],'b','LineWidth',1)
  // LI
  plot([Xa(7) Xa(11)], [Ya(7) Ya(11)],'b','LineWidth',1)

  // 4th Retangle - Microstrip line 
  // LS
  plot([Xa(6) Xa(7)], [Ya(6) Ya(7)],'b','LineWidth',1)
  //hold on
  // LE
  plot([Xa(6) Xa(5)], [Ya(6) Ya(5)],'b','LineWidth',1)
  // LD
  plot([Xa(7) Xa(8)], [Ya(7) Ya(8)],'b','LineWidth',1)
  // LI
  plot([Xa(5) Xa(5)], [Ya(8) Ya(8)],'b','LineWidth',1)

  //Clear feed lines
  // a
  dd=0.01; 
  plot([Xa(9)+dd Xa(6)-dd], [Ya(9) Ya(6)],'w','LineWidth',1)
  // b
  plot([Xa(6)+dd Xa(7)-dd], [Ya(6) Ya(7)],'w','LineWidth',1)
  // c
  plot([Xa(7)+dd Xa(11)-dd], [Ya(7) Ya(11)],'w','LineWidth',1)

  //Substrate
  plot([0 Aa], [0 0],'k','LineWidth',1)
  plot([0 Aa], [Aa Aa],'k','LineWidth',1)
  plot([0 0], [0 Aa],'k','LineWidth',1)
  plot([Aa Aa], [0 Aa],'k','LineWidth',1)
  diary off;