#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <time.h>
#include "mensagem.h"
#include "memvirtu.h"
#include "lowparam.h"
#include "lowmacro.h"
#include "lowsolid.h"
#include "eulerops.h"
#include "vectorop.h"
#include "genfunc_.h"
#include "analise_.h"
#include "shpshape.h"
#include "prop.h"
#include "mancommd.h"

real MSD_highNamePropriedadeArea(const char *name);
real MSD_highPropriedadeArea(Id sn);
double MSD_lowPropriedadeArea(SPTYPE s);
double MSD_lowPropriedadeAreaFace(FPTYPE f);
double MSD_lowPropriedadeAreaLaco(LPTYPE l);
real MSD_highNamePropriedadeVolume(const char *name);
real MSD_highPropriedadeVolume(Id sn);
double MSD_lowPropriedadeVolume(SPTYPE s);
double MSD_lowPropriedadeVolumeFace(FPTYPE f);
double MSD_lowPropriedadeVolumeLaco(LPTYPE l);
float MSD_SA(void);

void MSD_execManipulatePropriedade(void)
{
    int ip ;
    char onam[30] ;
    float area, volume, x1, x2, y1, y2, z1, z2, th, d1, d2, d3, d4, d5, d6, d7;
                                                   // definição das variáveis locais

    for (ip = 0 ; ip == 0 ; ) {
        switch(optin()) {                 // recupera as oções do comando

            // - Comando calcula a area - recupera o nome do solido
            case 'a':
                if (1 == sscanf(restbuf, "%s", onam)) {
                    area = MSD_highNamePropriedadeArea(onam) ;
                    ip = 1 ;
                    printf("area 10.03.2020 = %f\n", - area);
                }
                break ;

            // - Comando calcula o volume - recupera o nome do solido
            case 'v':
                if (1 == sscanf(restbuf, "%s", onam)) {
                    volume = MSD_highNamePropriedadeVolume(onam) ;
                    ip = 1 ;
                    printf("volume = %f\n", volume);
                }
                break ;

             // - Comando que calcula a rotação em relação a um eixo
            case 'r':
                if (8 == sscanf(restbuf, "%s %f %f %f %f %f %f %f", 
                                           onam, &x1, &y1, &z1, &z2, &y2, &z2, &th)) {
                //    MSD_highNamePropriedadeRotaciona(onam, x1, y1, z1, x2, y2, z2, th) ;
                    ip = 1 ;
                }
                break ;

            // - Comando para criar o primitivo barra L
            case 'l':
                if (6 == sscanf(restbuf, "%s %f %f %f %f %f", 
                                         onam, &d1, &d2, &d3, &d4, &d5)) {
                //    MSD_ highCreatePrimitivoL(onam, d1, d2, d3, d4, d5) ;
                    ip = 1 ;
                }
                break ;

            // - Outros comandos

        }
        if (ip == 0) {
            printf("-avrs nome do sólido\n") ;
            if (!lineins("? "))
                return ;
        }
    }
}

real MSD_highNamePropriedadeArea(const char *name)
{
    int sn ;                       // identificador do sólido

    if ((sn = MSD_getSolidIdFromName(name)) == -1) {
        fprintf(stderr, "area: nao encontrou o sólido %s!\n", name) ;
        return(ERROR) ;
    }
    return(MSD_highArea(sn)) ;
}

real MSD_highPropriedadeArea(Id sn)
{
    SPTYPE s ;                     // ponteiro para o sólido

    if ((s = MSD_getSolid(sn)) == SNIL) {
        fprintf(stderr, "area: nao encontrou o sólido %d!\n", sn) ;
        return(0.0) ;
    }
    return(MSD_lowPropriedadeArea(s)) ;
}

double MSD_lowPropriedadeArea(SPTYPE s)
{
    DPTYPE d;
    FPTYPE f;
    double area = 0;

    for (d = SolSShells(s) ; d != DNIL ; d = SheNextD(d))
        for (f = SheSFaces(d) ; f != FNIL ; f = FacNextF(f))
            area += MSD_lowPropriedadeAreaFace(f);
            MSD_lowSetNormal(s, TRUE) ;
MSD_lowSetEdgeAngle(s) ;
MSD_lowSetInfo(s) ;
    return(area);
}

double MSD_lowPropriedadeAreaFace(FPTYPE f)
{
    LPTYPE l;
    double area = 0;

    area = MSD_lowPropriedadeAreaLaco(FacFLOut(f));
    for (l = FacFLoops(f) ; l != LNIL ; l = LooNextL(l))
        if (l != FacFLOut(f))
            area -= MSD_lowPropriedadeAreaLaco(l);
    return(area);
}

/*** Calculo da area do laco ***/
double MSD_lowPropriedadeAreaLaco(LPTYPE l)
{
    HPTYPE he ;
    VPTYPE v1 ;
    vector aa, bb, cc, dd, vv1 ;

    veczer(dd) ;
    he = LooLEdg(l) ;
    v1 = HalVtx(he) ;
    he = HalNxt(he) ;
    do {
        veccopy(vv1, VerVCoord(v1)) ;
        vecminus(aa, VerVCoord(HalVtx(he)), vv1) ;
        vecminus(bb, VerVCoord(HalVtx(HalNxt(he))), vv1) ;
        cross(cc, aa, bb) ;
        vecplus(dd, dd, cc) ;
    } while ((he = HalNxt(he)) != LooLEdg(l)) ;
    return(-0.5 * dot(FacFeq(LooLFace(l)), dd)) ;
}

real MSD_highNamePropriedadeVolume(const char *name)
{
    int sn ;                       // identificador do sólido

    if ((sn = MSD_getSolidIdFromName(name)) == -1) {
        fprintf(stderr, "volume: nao encontrou o sólido %s!\n", name) ;
        return(ERROR) ;
    }
    return(MSD_highVolume(sn)) ;
}

real MSD_highPropriedadeVolume(Id sn)
{
    SPTYPE s ;                     // ponteiro para o sólido

    if ((s = MSD_getSolid(sn)) == SNIL) {
        fprintf(stderr, "volume: nao encontrou o sólido %d!\n", sn) ;
        return(0.0) ;
    }
    return(MSD_lowPropriedadeVolume(s)) ;
}

double MSD_lowPropriedadeVolume(SPTYPE s)
{
    DPTYPE d;
    FPTYPE f;
    double volume = 0;

    for (d = SolSShells(s) ; d != DNIL ; d = SheNextD(d))
        for (f = SheSFaces(d) ; f != FNIL ; f = FacNextF(f))
            volume += MSD_lowPropriedadeVolumeFace(f);
    
	MSD_lowSetNormal(s, TRUE) ;
    MSD_lowSetEdgeAngle(s) ;
    MSD_lowSetInfo(s) ;
    return(volume);
}

double MSD_lowPropriedadeVolumeFace(FPTYPE f)
{
    LPTYPE l;
    double volume = 0;

/*    area = MSD_lowPropriedadeAreaLaco(FacFLOut(f)); */
    for (l = FacFLoops(f) ; l != LNIL ; l = LooNextL(l))
/*        if (l != FacFLOut(f)) */
            volume += MSD_lowPropriedadeVolumeLaco(l);
    return(volume);
}

/*** Calculo da area do laco ***/
double MSD_lowPropriedadeVolumeLaco(LPTYPE l)
{
    HPTYPE he ;
    VPTYPE v1 ;
    vector aa, bb, cc, vv1 ;
    double dd = 0;

/*    veczer(dd) ; */
    he = LooLEdg(l) ;
    v1 = HalVtx(he) ;
    he = HalNxt(he) ;
    do {
        veccopy(vv1, VerVCoord(v1)) ;
        veccopy(aa, VerVCoord(HalVtx(he))) ;
        veccopy(bb, VerVCoord(HalVtx(HalNxt(he)))) ;
        cross(cc, aa, bb) ;
        dd += dot(vv1, cc) ;
    } while ((he = HalNxt(he)) != LooLEdg(l)) ;
    return(dd / 6.0) ;
}

double f(double y)
{
  if (y > -0.5 && y < 0.5)
    return y*y - 1000.0;
  return y * y + 5.0 * y + 6.0;
}

float MSD_SA(void)
{
    double s0 = -5.0;
    int n = 1;       // number of variables
    int p = n * 20;  // máximo de iterações local (perturbações da solução atual)
    int l = n * 10;  // Número máximo de sucessos por iteração
    double alfa = 0.98;
    double T0 = 5.0;
    double deltay = 5.0;

    double s = s0;
    double T = T0;
    int j = 1; // contador global
    int k = 1; // contador de soluções aceitas
    int c = 1; // crystallization factor
    srand(time(NULL));

    // loop principal
    while (k != 0) {
        int i = 1; // contador local
        k = 0;
        while (k < l && i < p) {
            double d = 0.0;
            for (int kk = 0 ; kk < c ; kk ++) {
                d += 2.0 * rand()/(double)RAND_MAX - 1.0;
            }
            d = d / (double)c * deltay;
            double delta = f(s + d) - f(s);

            if (delta < 0 || exp(-delta/T) > rand()/(double)RAND_MAX) {
                s += d;
                k ++;
                if (c > 1) c --;
            }
            else
                if (c < 80) c ++;
            i ++;
        }
        
        T *= alfa;
    }

    return s;
}